Analysis of anisotropic flow with Lee- Yang zeroes 
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We present a new method to extract anisotropic flow in heavy ion collisions from the genuine 
correlation among a large number of particles. Anisotropic flow is obtained from the zeroes in the 
complex plane of a generating function of azimuthal correlations, in close analogy with the theory 
of phase transitions by Lee and Yang. Flow is first estimated globally, i.e., averaged over the phase 
space covered by the detector, and then differentially, as a function of transverse momentum and 
rapidity for identified particles. The corresponding estimates are less biased by nonflow correlations 
than with any other method. The practical implementation of the method is rather straightforward. 
Furthermore, it automatically takes into account most corrections due to azimuthal anisotropics in 
the detector acceptance. The main limitation of the method is statistical errors, which can be 
significantly larger than with the "standard" method of flow analysis if the flow and/or the event 
multiplicities are too small. In practice, we expect this to be the most accurate method to analyze 
directed and elliptic flow in fixed-target heavy-ion collisions between 100 MeV and 10 GeV per 
nucleon (at the Darmstadt SIS synchrotron and the Brookhaven Alternating Gradient Synchrotron), 
and elliptic flow at ultrarelativistic energies (at the Brookhaven Relativistic Heavy Ion Collider, and 
the forthcoming Large Hadron Collider at CERN). 

PACS numbers: 25.75.Ld, 25.75.Gz, 05.70.Fh 



I. INTRODUCTION 

Study of anisotropic flow of particles [1] produced in 
relativistic heavy-ion collisions has emerged as an impor- 
tant tool to probe the early history, especially the ther- 
malization, of the dense fireball produced in these col- 
lisions. Anisotropic flow means that the azimuthal dis- 
tribution of particles produced in non-central collisions, 
measured with respect to the direction of impact param- 
eter, is not flat. This is characterized by the Fourier 
harmonics [2]: 

v n = (cosn((j) - $ R )) , (1) 

where <fi denotes the azimuthal angle of an outgoing par- 
ticle, is the azimuthal angle of the impact parameter 
is also called the orientation of the reaction plane), n 
is a positive integer, and angular brackets denote an av- 
erage over many particles belonging to some phase-space 
region, and over many events. Throughout this paper, 
we assume that colliding nuclei are spherical and parity is 
conserved: then, the system is symmetric with respect to 
the reaction plane and (sinn(0 — &r)) vanishes for all n. 
The two lowest harmonics v% and vi are named directed 
and elliptic flows. They have been studied extensively 
over the last several years, with reference to the wealth 
of data produced at GANIL [3], at the Darmstadt SIS 
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synchrotron [4-7] , the Brookhaven Alternating Gradient 
Synchrotron (AGS) [8-14], and the Super Proton Syn- 
chrotron (SPS) at CERN [15-19], as well as the new and 
upcoming data from the Relativistic Heavy Ion Collider 
(RHIC) at Brookhaven [20-25]. 

Since the reference direction is not known experi- 
mentally, measuring v n is a difficult task. An alternate 
idea is to extract it from the interparticle correlations 
which arise indirectly due to the correlation of each par- 
ticle with the reaction plane. However, standard methods 
of correlating a particle with an estimate of the reaction 
plane [26-28], or of correlating two particles with each 
other [29] were shown to be inadequate at ultrarelativis- 
tic SPS energies [30] due to the smallness of the flow, and 
the comparatively large magnitude of "nonflow" correla- 
tions [28, 31] due to transverse momentum conservation, 
resonance decays, etc., which are usually neglected in 
these methods. At the higher RHIC energies, it was ar- 
gued that nonflow correlations due to minijets could even 
dominate the measured correlations [32]. 

These shortcomings of conventional methods moti- 
vated the development of new methods based on a cu- 
mulant expansion of multiparticle correlations [33-35]. 
The cumulant of the fc-particle correlation, where k is a 
positive integer, isolates the genuine /c-particle correla- 
tion by subtracting the contribution of lower-order cor- 
relations. Nonflow correlations, which generally involve 
only a small number of particles (typically, two or three 
in the case of resonance decays, possibly more in the case 
of minijets) contribute little to the cumulant if k is large 
enough. Note that this is also true for correlations from 
global momentum conservation, although the latter in- 
volves all particles [36]. On the other hand, anisotropic 
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flow is a genuine collective effect, in the following specific 
sense: in a given event, all azimuthal angles are corre- 
lated with the reaction plane azimuth which varies 
randomly from one event to the other. In the laboratory 
frame, where is unknown, this results in azimuthal 
correlations between all particles, or at least a significant 
fraction of the particles. Therefore, unlike nonflow cor- 
relations, anisotropic flow contributes to cumulants of all 
orders k. Neglecting the contribution of nonflow correla- 
tions to the cumulant, one thus obtains an estimate of v n , 
which was denoted by v n {k} in Ref. [34], and becomes 
more and more reliable as k increases. 

Cumulants of four-particle correlations were first used 
to measure elliptic flow at RHIC [24]. However, it was 
argued that experimental results could still be explained 
by nonflow correlations alone [37] at this order. Higher- 
order cumulants, of up to 8 particles, were constructed 
at SPS [18], and provide the first quantitative evidence 
for anisotropic flow at ultrarelativistic energies. In prac- 
tice, the cumulant of 8-particle correlations is obtained 
by evaluating numerically the 8 th derivative of a gener- 
ating function. This is rather tedious and numerically 
hazardous. 

In this paper, we propose to study directly the large- 
order behavior of the cumulant expansion, rather than 
computing explicitly cumulants at a given order. Corre- 
lating a large number of particles is the most natural way 
of studying genuine collective motion in the system. Fur- 
thermore, finding the large-order behavior turns out to be 
simpler in practice than working at a given, finite order. 
As we shall see, the large-order behavior is determined by 
the location of the zeroes of a generating function in the 
complex plane [38]. In this respect, our method is in close 
analogy with the theory of phase transitions formulated 
50 years ago by Lee and Yang [39, 40]: at the critical 
point, long-range correlations appear in the system; as 
a consequence, the zeroes of the grand partition func- 
tion come closer and closer to the real axis as the size of 
the system increases. A similar phenomenon occurs here, 
and anisotropic flow appears as formally equivalent to a 
first order phase transition. 

The idea of studying anisotropic flow through the 
correlation of a large number of particles is not new. 
The same idea underlies global event analyses through 
a three-dimensional sphericity tensor [41], which led to 
the first observation of collective flow at Bevalac [42]. 
A similar two-dimensional analysis, restricted to the 
transverse plane [2, 43], led to the discovery of flow 
at the Brookhaven Alternating Gradient Synchrotron 
(AGS) [9]. With these methods, however, one could not 
study anisotropic flow differentially, that func- 
tion of transverse momentum and rapidity for identified 
particles. This is why they were soon superseded by the 
more detailed, although less reliable, two-particle meth- 
ods. One of our important results is that the method pre- 
sented in this paper also allows one to analyze differential 
flow, and its implementation is rather straightforward. 

The paper is organized as follows. In Sec. II we 



present in a self-contained manner our recipes to calcu- 
late the "asymptotic" integrated and differential flows. 
The derivations of these results are given in Sec. Ill, 
where they are related with the Lee- Yang theory. The 
remaining Sections are devoted to detailed discussions of 
the various sources of error. It is now well known that 
(even when anisotropic flow is absent) "nonflow" corre- 
lations may give a spurious result when the various tech- 
niques of flow analysis are applied. In Sec. IV, we discuss 
the magnitude of this spurious flow, and show that it is 
significantly smaller than with any previous method: the 
main point of this paper is that the present method is the 
most efficient way to disentangle anisotropic flow from 
other effects. When anisotropic flow is present in the 
system, nonflow correlations produce a systematic error, 
which is estimated in Sec. V. The effect of fluctuations of 
flow (due for instance to variations of impact parameter 
in a given centrality bin) on our results is discussed in 
Sec. VI. Statistical uncertainties on the flow estimates 
within the method are computed in Sec. VII. They are 
the main practical limitation of our method. Acceptance 
corrections due to limited azimuthal coverage of the de- 
tectors are discussed in Sec. VIII. Section IX is a sum- 
mary. Finally, four appendices are devoted to further 
discussions and calculations. 



II. RECIPES FOR EXTRACTING GENUINE 
COLLECTIVE FLOW 

In this Section, we show how to obtain an estimate of 
the flow v n from the genuine correlation among a large 
number of particles. The fact that azimuthal angles of 
outgoing particles are correlated with the azimuth of the 
reaction plane allows one to construct in each event 
a reference direction, which is called an "estimate of 
the reaction plane" in the standard Danielewicz-Odynicc 
method [26]. This direction is defined by the flow vector 
of the event, as recalled in Sec. II A. Note that, unlike 
other multiparticlc methods [34] , the present one is based 
on the flow vector. 1 

The flow analysis then proceeds in two successive steps. 
The first one is to estimate how the flow vector is cor- 
related with the true reaction plane. More precisely, 
we estimate the mean projection of the flow vector on 
the true reaction plane. This quantity is a weighted 
sum of the individual flows v n of all particles over phase 
space, which we call "integrated flow" and denote by V n . 
The method used to estimate V n experimentally is de- 
scribed in Sec. II B. This first step is the equivalent of 
the subevent method in the standard flow analysis, from 
which one estimates the event-plane resolution. Here 



1 In Appendix A, we propose an alternate version of the present 
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however, as well as in the cumulant method based on 
the flow vector [33], subevents are not needed. 

The second step in the analysis is to use this reference 
integrated flow to analyze "differential flow," i.e., flow in 
a restricted phase-space window (e.g., as a function of 
transverse momentum and rapidity for a given particle 
type), which is the goal of the flow analysis. When an- 
alyzing a p-th harmonic of the differential flow (v p ), the 
flow vector can be chosen in the same harmonic p, or in a 
lower harmonic n whose p is a multiple. One chooses the 
option which yields the best accuracy. 2 For instance, 
differential elliptic flow v-i can be analyzed using as a 
reference either the integrated elliptic flow (at ultrarel- 
ativistic energies [15, 17-24]) or the integrated directed 
flow (at lower energies [5, 14, 44]). 

The method to extract differential flow is described in 
Sec. IIC. Following the notations of Refs. [33, 34], differ- 
ential flow in the Fourier harmonic p will be denoted by 
v' p . The estimates of V n and v' p obtained with the present 
method will be denoted by l^ijoo} and v' p {co}, respec- 
tively, where the symbol oo means that it corresponds to 
the large-order behavior of the cumulant expansion, as 
will be shown in Sec. III. 

Section II D discusses briefly acceptance issues arising 
when the detector does not have full azimuthal coverage. 
Section HE shortly deals with statistical errors, which 
are the main limitation of our method. Both issues are 
discussed in more detail in Sec. VIII and Sec. VII, re- 
spectively. 

A. The flow vector 

The first step of the flow analysis is to evaluate, for 
each event, the flow vector of the event. It is a two- 
dimensional vector Q = (Q x ,Q y ) defined as 

M 

Qx = ^2 w j cos(n^) 
i=i 

M 

Qy = ^2 w j sin(n^), (2) 
i=i 

where n is the Fourier harmonic under study (n = 1 for 
directed flow, n — 2 for elliptic flow), the sum runs over 
all detected particles, M is the observed multiplicity of 
the event, (f>j are the azimuthal angles of the particles 
measured with respect to a fixed direction in the labora- 
tory. 

The coefficients Wj in Eq. (2) are weights depending on 
transverse momentum, particle mass and rapidity. The 



2 Note also that the only way to determine the sign of v p , with 
p > 2, is to use a flow vector in a lower harmonic n < p. Thus, at 
SPS the sign of elliptic flow V2 is determined using the reference 
from directed flow vi, while its magnitude is determined more 
accurately using the flow vector in the second harmonic [18] . 



best choice of Wj is that leading to the smallest statistical 
error, by maximizing the flow signal. The weight should 
ideally be proportional to the flow itself, Wj(pT,y) oc 
v n(PT,y) [33]. Otherwise, reasonable choices are w oc 
y — ycm, i-c, the rapidity in the center-of-mass frame, for 
directed flow, and w — pt for elliptic flow. 

When comparing events with different multiplicities 
M, one may in addition have weights depending on M. 
Weights proportional to \j\f~M were used in Ref. [33], 
and proportional to 1/M in Ref. [34], to minimize the 
effect of multiplicity fluctuations. This complication is 
unnecessary here. This issue is discussed in Appendix B. 

The discussions in this paper will be illustrated by ex- 
plicit numerical and analytical examples, in which wc 
choose unit weights Wj = 1, for the sake of simplicity. 

The flow vector was first introduced in Ref. [26] . The 
azimuthal angle n$„ of Q is conventionally used in the 
analysis of anisotropic flow [28] in order to estimate the 
orientation of the reaction plane of the event. 

Our method uses the projection of Q on a fixed, arbi- 
trary direction making an angle nO with respect to the 
x-axis. We denote this projection by Q e : 

Q 6 = Qx cos(n6) + Q y sin(n#) 

M 

= *^2 W 3 cos(n(<^- - 9)). (3) 
j=i 

As we shall see, the whole flow analysis can in principle 
be performed using this projection of the flow vector on 
a fixed direction. One thus obtains an estimate of the 
integrated flow V®{oo} (Sec. II B), which is then used as 
a reference to derive an estimate of the differential flow 
*0°o} (Sec. IIC). 

In practice, however, one should perform the anal- 
ysis for several equally spaced values of 9, typically 
9 = (k/p)(n/n) with k = 0, . . . ,p— 1 andp = 4 or 5. This 
gives several values of V^{oo} and w^ n {oo}, which are 
then averaged over 9. This yields our final estimates of 
integrated flow, V n {oo} 1 and differential flow, v' mn {oo}. 
As we shall see in Sec. VII, they have smaller statistical 
error bars than each individual V^f {oo} and w^ n {oo}. 

B. Integrated flow 

Integrated flow is defined as the average value of the 
flow vector projected on the unit vector with angle n&R 
(for n = 1, this is the reaction plane): 

V n = (Q x cos(n$ K ) + Q y sin(n$ fl )) 

= (4) 

where we have used the notation introduced in Eq. (3). 
This quantity was denoted by Q in Ref. [27], and by (Q) 
in Ref. [33]. Note that unlike the flow coefficients v ni 
which are dimensionless, our integrated flow involves the 
weights present in the flow vector (2), which can have 
a dimension. Since the flow vector definition involves a 
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sum over all particles, the integrated flow scales like the 
multiplicity M. With unit weights, 



V„ = Mv n , 



(5) 



where v n in the right-hand side (rhs) is to be understood 
as an average over the phase space covered by the detec- 
tor acceptance, and we have neglected fluctuations of the 
multiplicity M for simplicity. 

Let us now explain how an estimate of V n can be ob- 
tained in a real experiment. For a given value of 9, we 
first define the following generating function, which de- 
pends on an arbitrary complex variable z: 



(6) 



where curly brackets denote an average over a large num- 
ber of events, N evts , with (approximately) the same cen- 
trality. This generating function has the symmetry prop- 
erties 



G e+7r/n (z) = G°(-z) 
(following from Q 0+7T / n = -Q 6 ), and 
[G e {z)\* = G e (z*), 



(7) 



(8) 



where the star denotes complex conjugation. The second 
identity simply expresses that G s (z) is real for real z. An 
alternative choice of the generating function is presented 
in Appendix A. 

In order to obtain the integrated flow, one must evalu- 
ate G e (z) for a large number of values of z on the upper 
half of the imaginary axis (that is, z — ir with r real 
and positive). One must then take the modulus [recall 
that G e (z) is a complex number], |G 9 (ir)|, and plot it 
as a function of r. On the imaginary axis, the symme- 
try properties, Eqs. (7) and (8) translate into the follow- 
ing relations: G e+7T / n (ir) = G e (-ir) = [G e (ir)}*. Tak- 
ing the modulus, this yields |G 9+7r/n (ir)| = \G e (~ir)\ = 
\G (ir)\. This identity shows that 9 and 9 + Ti/n are 
equivalent, so that one can restrict 9 to the interval 
[0,7r/ra]. It also shows that z = ir and z — —ir are 
equivalent, which is the reason why we restrict ourselves 
to positive r values. 

To illustrate the variation of | G e (ir) | as a function of 
r, we have computed it for simulated data. The data set 
contained iVevts = 20000 events with multiplicity M — 
300. Each event consisted of 10 "differential bins" with 
equally spaced elliptic flow values ranging from 4.2 to 
7.8%, resulting in an average «2 = 6%, and with a higher 
flow harmonic i>4 = 3% for all particles. The flow vector, 
Eq. (2), was constructed using unit weights Wj = 1 and 
n = 2. We assumed that the detector acceptance had 
perfect azimuthal symmetry. The function |G e (ir)| is 
shown in Fig. 1 for 9 = 0. Due to rotational symmetry 
(for a perfect detector), the behavior would be similar 
for another value of 9, up to statistical fluctuations. The 
function |G e (ir)| starts at a value of 1 for r = 0, and it 
then decreases and oscillates. 
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FIG. 1: Variation of \G e (ir)\ with r, for iV cvt s = 20000 sim- 
ulated events, with M = 300 particles per event and a mean 
elliptic flow V2 = 6%, to match the typical value in a mid- 
central Au+Au collision at ^/s^7 = 130 GeV, as analyzed 
by the STAR Collaboration at RHIC [23]. The crosses are 
the values of |G (ir)\ for 8 = 0. The solid line displays the 
expected value |G c .i.(ir)| [defined by Eq. (23)]. 



Our estimate of integrated flow is directly related to 
the first minimum of |G 9 (ir)|. Let us denote by rg the 
value of r corresponding to the first minimum. The cor- 
responding estimate of integrated flow is defined as 



i,«r i _ Joi 



(9) 



where joi — 2.405 is the first root of the Bessel function 
Jo(x). This result will be justified in Sec. III. As one can 
see in Fig. 1, the variation of |G e (ir)| near its minimum 
is quite steep. Therefore, from the numerical point of 
view, one should rather determine the minimum of the 
square modulus | G e (ir) | . 3 

As will be shown in Sec. Ill, the value of |G e (ir)| at its 
minimum would be zero in the limit of an infinite number 



The following procedure can be used. One first computes 
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where V ma x is some a priori upper bound on the expected value 
of the integrated flow V n , e is a small increment, and N is an 
integer varying between and V mm /e — 1. Next, denote by ./Vo 
the first value of N for which > f%, so that /jy < f^ fQ _ 1 

and fff a < Then, the integrated flow is approximately 

given by 
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of events. One can check experimentally that it is com- 
patible with zero within expected statistical fluctuations. 
Anticipating on Sec. VII C [see the discussion following 
Eq. (72], the following inequality should hold with 95% 
confidence level: 



|G e K)| < 



(10) 



where -/V cvt s is the number of events used in the analysis. 

Eventually, the final estimate V„{oo} is obtained by 
averaging V® {oo} over 9: 



F4oc}^-y>„ fc ^"{oc}, 



p 



(11) 



fc=0 



with p typically 4 or 5. 

This procedure was applied to the simulated data used 
for Fig. 1. As stated above, we used constant unit weights 
u>j = 1 in Eq. (2), so that Eq. (5) holds. Estimates 
V 2 e {oo} were then computed for many values of 9 (for 
the sake of illustration, we did not restrict ourselves to 4 
or 5 values, as would have been enough, see Sec. VII E). 
Since we assumed that the detector is perfect, V® {°°} 
should be independent of 9, up to statistical fluctua- 
tions. The results of the analysis are displayed in Fig. 2, 
where V® {oo}/M is plotted as a function of 9. The 9- 
dependence is smooth and has a small amplitude. The 
estimates for 9 = and 9 = ir/2 coincide, as expected 
from the (7r/n)-periodicity. 
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FIG. 2: Reconstruction of integrated flow. The simulated 
data and the definition of the flow vector are the same as in 
Fig. 1. The plot shows V${oo}/M as a function of 9. The 
lines display the expected value with statistical error bars 
calculated as explained in Sec. VII E. 



The salient feature in Fig. 2 is the accuracy of the re- 
construction of the input. As will be shown in Sec. VII E, 




the expected statistical uncertainty on each individual 
estimate V® {oo}/M is 0.092%. One sees in the figure 
that most values of V% f{oo}/M are within error bars. 
The expected statistical error bar on the average over 9, 
V 2 {oo}/M, is 0.051% (see Sec. VII E), i.e., smaller than 
the statistical error on each V/joo} by almost a factor of 
2. We obtain V2{oo}/M = 5.95%, in perfect agreement 
with the input value v-i = 6%. 

Finally, let us mention that, like every other method of 
flow analysis, the present one cannot determine the sign 
of integrated flow: the estimate given by Eq. (9), and 
its average over 9, are positive by definition, while the 
integrated flow defined by Eq. (4) can be negative. The 
reason is that our procedure is unchanged if the sign of 
the flow vector is changed in all events [this amounts to 
changing z into — z in Eq. (6)], while this transformation 
changes the sign of the integrated flow V n . Strictly speak- 
ing, V^{oo} should therefore be considered an estimate 
of |V^|, rather than V n . The sign must be determined 
independently, or assumed. 



C. Differential flow 

Using an estimate of integrated flow in the Fourier har- 
monic n, one can analyze differential flow in harmonics 
which are multiples of n, i.e., mn, where m is an arbi- 
trary integer. Following the notations of Ref. [34] , we de- 
note by v' mn the differential flow corresponding to a given 
phase-space window in this harmonic. We call "proton" 
any particle belonging to the phase-space window under 
study, and denote by ip its azimuthal angle. 

For a given value of the angle 9, the corresponding 
estimate of differential flow v' mn is given by: 



t {oo} 



Ji(j 01 ) R / {cos[mn(^ - 6»)] e^"} 



i™- 1 {Q e e<Q e } 



(12) 

where Re denotes the real part, and r® has been defined 
in Sec. II B. Two different sample averages appear in the 
rhs of this equation: the numerator is an average over 
all "protons" in all events, while the denominator is an 

e e ir $ Q B 



average over events. Note that the term ^Q e e lT ' e ° ( ^ e j in 

the denominator is in fact the derivative of G (z) with 
respect to z [see Eq. (6)], evaluated at z = ir®. 

The numerical coefficient Ji(joi)/^m(joi) m Eq. (12) 
involves the ratio of two Bessel functions. It takes the 
values 1 for m — 1 and jai/2 ~ 1.202 for m — 2. In the 
case m = 1 (lowest harmonic) , one recovers the estimate 
of integrated flow V®{oo} by integrating the correspond- 
ing estimate of differential flow v'^{oo}, over all phase 
space, which amounts to summing over ip in the rhs of 
Eq. (12), with the appropriate weighting. 

If cos is replaced by sin in the numerator of the rhs 
of Eq. (12), the result should be zero within statistical 
errors if parity is conserved. This can easily be checked 
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experimentally. A non-zero result could be a signature of 
parity violation [45, 46]. This issue will not be discussed 
further in this paper. 

As in the case of integrated flow, the estimates of differ- 
ential flow given by Eq. (12) have a periodicity property, 

namely ^"''"{oo} = u^ n {oo}. Indeed, changing 6 into 
+ n/n amounts to replacing the term in brackets by its 
complex conjugate, which does not affect the value of the 
real part. As above, the estimates u^ n {oo} must be av- 
eraged over 9 in the interval [0, ir/n], in order to obtain 
an estimate with a reduced statistical uncertainty. 

Please note that there appear "autocorrelations" in the 
numerator of Eq. (12): one correlates the angle ?/> with 
Q e , which itself involves in general the angle ip [since 
the summation in Eq. (2) runs over all detected parti- 
cles]. In the standard method of flow analysis [26, 28], 
autocorrelations are large, so that one has to remove the 
particle with angle ip from the flow vector. Here, how- 
ever, autocorrelations do not produce a spurious flow by 
themselves, as will be explained in Sec. HIE 4. For the 
lowest harmonic v' n , they lead to a very small correc- 
tion and need not be subtracted. 4 For higher harmonics 
v' 2n , v' 3n ..., errors due to autocorrelations are larger so 
that one may prefer to subtract them. However, errors 
of the same order of magnitude are generally expected 
from nonflow correlations. We come back to this issue in 
Sec. VB. 

Figure 3 displays the result of the analysis of the same 
Monte-Carlo simulation as in Figs. 1 and 2, using a de- 
tector with perfect azimuthal symmetry. We present the 
differential flow results in "bin 8," corresponding to input 
values of v 2 = 7%, v' 4 = 3%, and an average proton num- 
ber of approximately 30 per event, i.e., a total number of 
protons N' ~ 6 x 10 5 . Two different analyses were per- 
formed. The first one followed the procedure presented 
in this Section. In a second analysis, we corrected for 
autocorrelations, subtracting the contribution of the pro- 
ton with angle ip from the flow vector before multiplying 
cos[mn(ip — 9)] by e tr °Q in the numerator of Eq. (12). 
One sees in Fig. 3 that the v' 2 result is essentially insensi- 
tive to autocorrelations. 5 We shall see in Sec. VII F that 
the expected statistical errors on v' 2 {oo} and ^{oo} are 
0.47% and 0.57%, respectively. Values mostly fall within 
the expected range around the input value, except for v\ 
when autocorrelations are not subtracted. 

After averaging over 9, one obtains v' 2 {oo} = 7.00 ± 



0.26% and v' 4 {oo} = 3.60 ± 0.29% with autocorrelations 
and v' 2 {oo} = 7.20 ± 0.26%, v' 4 {oo} = 3.03 ± 0.29% 
when autocorrelations are removed, where statistical un- 
certainties have been computed with the help of formulas 
given in Sec. VII F. For the lower harmonic v' 2 , results are 
in good agreement with the input value v' 2 = 7%, whether 
or not autocorrelations are subtracted. For the higher 
harmonic v' Al a discrepancy with the input value v\ = 3% 
appears when autocorrelations are not subtracted. This 
error will be evaluated analytically in Sec. VB. 
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FIG. 3: Reconstruction of differential flow. The simulated 
data and the definition of the flow vector are the same as in 
Figs. 1 and 2. The plot shows v'-2{oo} and t>4 S {oo} as a func- 
tion of 8. Squares: v' 2 {oo} with autocorrelations. Diamonds: 
v'2 {oo} with autocorrelations removed. Crosses: V4{oo} with 
autocorrelations. Circles: v'floo} with autocorrelations re- 
moved. The solid lines display the expected value with sta- 
tistical error bars calculated as explained in Sec. VII F. 



Like other methods of flow analysis, the present proce- 
dure has a global sign ambiguity, due to the fact that the 
sign of the integrated flow V n cannot be reconstructed. 
In Eq. (12), both V^{oo} and 7q are positive. If the 
true integrated flow V n is also positive, then our estimate 
v'J^ n {oo} has the correct sign. If, on the other hand, V n is 
negative, then w^ n {oo} should be multiplied by (— l) m . 
[equivalently, one can change the sign of V®{oo} and 7q 
in Eq. (12)]. 



4 If autocorrelations are subtracted, in particular, one no longer 
recovers exactly the integrated flow by integrating the differential 
flow over phase space. 

5 A detailed study of systematic errors, carried out in Appendix 
D, explains why the flow values with autocorrelations subtracted 
[Eq. (DI5)] are slightly larger (instead of much smaller in the 
standard analysis) than when autocorrelations are not taken into 
account [Eq. (Dll)]. For the present simulation, e defined in 
Eq. (D5) equals approximately 0.012, which explains the relative 
difference of 2.5% between diamonds and squares in Fig. 3. 



D. Acceptance corrections 

The standard event-plane analysis requires that the 
flow vector Q has a perfectly isotropic distribution in 
azimuth. Since real detectors are not perfect, this re- 
quires one to use various flattening procedures [28]. One 
of the nice features of the cumulant expansion is that it 
isolates physical correlations by subtracting out the con- 
tribution of detector asymmetries [33]. The same occurs 
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here, so that flattening procedures are not needed. De- 
tector asymmetries can never produce a spurious flow by 
themselves with the method described above. This holds 
even if the detector has a very limited azimuthal cover- 
age. 

In Sec. VIII, we demonstrate that the effects of strong 
azimuthal asymmetries in a detector are twofold. First, 
the estimates V 2 e {oo} depend on 9 [see Eq. (105)]. Sec- 
ond, the average estimate T^{oo} does not exactly co- 
incide with the true value, but differs by a multiplica- 
tive coefficient (which also controls the ^-dependence of 
V 2 {oc}). In most cases, however, this coefficient is so 
close to unity that correcting for it is not even necessary. 

This is illustrated by Fig. 4: we analyzed a simu- 
lated data set with similar statistics as that of Figs. 1-3, 
namely 20000 events of 300 emitted particles, with an av- 
erage elliptic flow v 2 — 6% (but now a vanishing V4). We 
assumed that the detector had a blind angle of 60 de- 
grees, i.e., that one sixth of the azimuthal coverage is 
missing. We can clearly see the oscillation of V 2 e {oo} 
as a function of 9. After averaging over 9, however, we 
find V 2 {oo}/ (M) = 6.004 ± 0.066%, where (M) denotes 
the mean multiplicity of detected particles, 6 in perfect 
agreement with the input value v 2 = 6%. Note that 
the statistical error bar is slightly larger than in Fig. 2 
due to the reduced multiplicity ((M) = 250 instead of 
300). With this particular detector configuration, the 
analytical calculation presented in Sec. VIII shows that 
the reconstructed V2{oo} should be divided by 1.0068 in 
order to take detector effects into account [see Eq. (106)], 
a very small correction indeed. 








ji/8 it/4 3it/8 it/2 



As a consequence of the blind angle, the number of detected 
particles is smaller than the number of emitted particles, which 
was 300 particles for all events in our simulation, and fluctuates 
around the mean value (M) = 250. 



FIG. 4: Reconstruction of integrated flow, with a detector 
which does not detect particles with 150° < 4> < 210°. The 
simulated data and the definition of the flow vector are the 
same as in Fig. 1. The plot shows V% '{00} / (M) as a func- 
tion of 8. The lines display the theoretical value taking into 
account detector effects, Eq. (105), with statistical error bars. 



Similarly, acceptance inefficiencies have the same two 
effects on the differential flow estimates obtained with the 
present method. Thus, i^„{oo} becomes ^-dependent, 
while v' mn {co} differs from v' mn by a multiplicative factor, 
see Sec. VIII. 

Finally, let us emphasize that since the method auto- 
matically takes into account most detector asymmetries, 
it is important to take as input in the construction of the 
generating function, Eq. (6), the azimuthal angles of the 
particles exactly as they are measured in the laboratory: 
flattening procedures are useless here, and might even 
bias the analysis. 



E. When is the method applicable? 

As we shall demonstrate in the following Sections, the 
main limitation of the method arises from statistical er- 
rors, while systematic errors inherent to the method are 
under good control. In Sec. VII, we shall see that statisti- 
cal errors depend on the so-called "resolution parameter" 
X- Referring the reader to Sec. VII A for further details 
on this parameter, and in particular on its experimen- 
tal determination, we simply want to recall here that it 
is roughly given by \ — v n \f~M, where M is the typi- 
cal number of detected particles per event, and v n the 
typical magnitude of the flow. 

Compared to other methods of flow analysis, statis- 
tical errors depend very strongly on x- This is because 
our method involves the correlation between a large num- 
ber of particles. It is therefore important to use as many 
particles as possible in constructing the flow vector, so as 
to maximize x- I n particular, since our method is very 
stable with respect to effects of detector acceptance as 
explained above, one should not do any cuts in phase 
space. Furthermore, since the method is also stable 
with respect to "nonflow" correlations, one should not 
worry about possible "double countings" due to multi- 
ple hits or showering effects in the detector. Most cur- 
rent heavy ion experiments consist of several detectors 
of various types (calorimeters, tracking chambers) cover- 
ing various regions of phase space (central rapidity, for- 
ward/backward rapidity). It is not common in the flow 
analysis to combine the information from different detec- 
tors in constructing the flow vector of the event. Here, it 
is important to do so: even a few additional particles may 
result in a significant decrease of the statistical error. Fi- 
nally, properly chosen weights can significantly increase 
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the value of the resolution parameter x- 7 

An over-simplified rule of thumb, which emphasizes the 
role of x would be the following: 

• if x > 1; the statistical error on the flow is not 
significantly larger than with the standard method 
(by a factor of at most 2, see Tables III, VI and 
VII), while systematic errors due to nonflow effects 
are much smaller. Thus, the present method should 
be used, and statistics will not be a problem; 

• if 0.5 < x < I, the method can be used, but it re- 
ally is most important to optimize weights, so as to 
increase x, possibly by performing two analyses of 
the same data sets: in the first analysis, adopting 
(educated) guesses for the weights; and in the sec- 
ond pass, using as weights Wj the differential flow 
values obtained in the first analysis; 

• if x < 0.5, statistical errors are too large, and 
the present method cannot be used; increasing the 
number of events barely helps here; in this case, one 
should use the cumulant methods of Refs. [34, 35] , 
which still apply if the number of events is large 
enough [18]. 



III. COLLECTIVE FLOW AND ZEROES OF 
THE GENERATING FUNCTION 

Let us now justify the procedure presented in the pre- 
vious Section. We begin with a detailed discussion of our 
hypotheses (Sec. Ill A). In Sec. Ill B, we define the cumu- 
lants of the distribution of Qf . We argue that large-order 
cumulants isolate collective effects from other, nonflow 
correlations, which generally involve only a few particles. 
We then show that large-order cumulants are uniquely 
determined by the location of the zeroes of G 6 (z) in the 
complex plane of the variable z. 

Anisotropic flow is a particular case of collective be- 
havior, which can be analyzed using this procedure. For 
this purpose, we compute the generating function G (z) 
when anisotropic flow is present (Sec. IIIC). This al- 
lows us to relate the position of the first zero of G e (z) 
to the integrated flow V n . The analogy of this approach 
with the Lee- Yang theory of phase transitions is shown 
in Sec. HID. 

The discussion of Sees. Ill B and III C is generalized 
to differential flow in Sec. HIE. The relation of this ap- 
proach with the cumulant expansion of Refs. [33, 34] is 
described in Appendix B. 



7 For instance, in the NA49 flow analysis, the values of the pt- 
wcightcd integrated elliptic flow are 20% larger than the corre- 
sponding values obtained with unit weights [18]. 



A. Preliminaries 

The derivation of our results relies on the following 
four hypotheses: 

• A) the event multiplicity, M, is much larger than 
unity; 

• B) for a fixed orientation of the reaction plane 
each particle in the system is correlated only to 
a small number of particles, which does not vary 
strongly with nuclear size and impact parameter; 

• C) the number of events available, V ovts , is arbi- 
trarily large; 

• D) the detector is azimuthally symmetric. 

Let us comment briefly on these hypotheses. With hy- 
pothesis C), all sample averages [such as {exp(zQ )} in 
Eq. (6)] become true statistical averages, which we write 
with angular brackets (i.e., (exp(zQ e ))). Hypothesis A) 
is verified in practice for heavy-ion collisions at high en- 
ergies. Hypothesis D) is not essential, but simplifies the 
calculations. Together with C), it implies for instance 
that the generating function, G e (z), is independent of 9. 

Hypothesis B) is the crucial assumption, where much 
physics is hidden, so let us comment on it in a little more 
detail. The idea is that for a fixed geometry (same impact 
parameter, same orientation in space of the colliding nu- 
clei), correlations are of the same order of magnitude as 
if the nucleus- nucleus collision was a superposition of in- 
dependent nucleon-nucleon collisions. This is a very rea- 
sonable assumption at ultrarelativistic energies 8 if the 
sample of events used in the analysis have exactly the 
same geometry, i.e., exactly the same impact parameter 
and <&r. Then, nucleon-nucleon collisions occurring in 
different places in the transverse plane are uncorrelated, 
simply by causality: the transverse size is much larger 
than the time scale of the collision due to Lorentz con- 
traction. Final state interactions may of course induce 
or (more likely) destroy correlations but we assume that 
they remain short-ranged, in the sense that they involve 
only a few particles, whose number remains roughly con- 
stant as the system size (i.e., the number of participating 
nucleons) increases. 9 This hypothesis breaks down in the 
vicinity of a phase transition, where critical fluctuations 
induce long-range correlations [47]. This would be even 
more interesting than anisotropic flow itself. 

Since the overall number of particles is large [hypoth- 
esis A)], hypothesis B) implies that each event can be 



We could not find a similar argument at lower energies, but nev- 
ertheless, we feel that the assumption remains valid. Anyway, 
one should be aware that it underlies all methods of collective 
flow analysis. 

Note that we do not exclude the possibility that there exist long- 
range correlations in phase space, for instance in rapidity. 
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split in some way into a large number of independent 
"subevents," whose number N roughly scales with the 
system size, i.e., with the multiplicity M. 10 In prac- 
tice, this means that if one creates an artificial event 
by taking the first subevent from one event, the second 
subevent from another event, etc., the resulting "mixed" 
event looks exactly like a real event (note that we have 
assumed that all events have exactly the same reaction 
plane therefore such mixed events cannot be con- 
structed experimentally: we use them simply as an image 
to illustrate our hypothesis). 

If N denotes the number of independent subevents, 
then Q can be written as the sum of N independent 
contributions: Q e = J2jLi Qj- We may write 

(e zQ6 \* R ) = f[(e* Q °\Z R ), (13) 

where the notation (. . . \<&r) denotes an average value 
taken over many events having exactly the same reaction 
plane Qr. The logarithm of this expression scales like 
the number of independent subsystems N, which itself 
scales like the multiplicity M: 

In (e zQ6 ~ M(az + bz 2 + cz 3 ■ ■ ■), (14) 

where the coefficients in the expansion are independent 
of the system size, and of order unity if the weights in 
Eq. (3) are of order unity. This equation is the mathe- 
matical formulation of hypothesis B), on which the fol- 
lowing discussion relies. Note that global momentum 
conservation, although it involves all particles, effectively 
behaves as a short-range correlation [36] , so that it does 
not invalidate Eq. (14). If the azimuthal distribution is 
symmetric (i.e., no anisotropic flow, azimuthally sym- 
metric detector), the left-hand side (lhs) is independent 
of 8 by azimuthal symmetry, and Eq. (7) then implies 
that it is an even function of z. As a consequence, odd 
terms vanish in the rhs of Eq. (14). 

B. Cumulants 

The cumulants (({Q e ) k )) are defined as the coefficients 
in the power series expansion of \nG e (z) [48]: 

+ 00 

fc=l K - 

If there is no anisotropic flow, outgoing particles are not 
correlated to the reaction plane Then, both sides 



of Eq. (14) are independent of The lhs is simply 

InG (z) according to the definition, Eq. (6). One con- 
cludes that the cumulants scale linearly with the size of 
the system, i.e., with the total multiplicity M. 

This scaling law breaks down if there is anisotropic 
flow, in which case the rhs of Eq. (14) depends on $>r. It 
also breaks down if hypothesis B) is not satisfied, i.e., if 
there are other collective effects (other than global mo- 
mentum conservation [36]) in the system. Then, cumu- 
lants generally scale with the total multiplicity like M k . 
This scaling law is most natural: Q 6 in Eq. (3) is the sum 
of M terms, so that in the generating function, Eq. (6), 
z is always multiplied by M terms. It is therefore nat- 
ural that z k in Eq. (15) goes with a factor proportional 
to M k . This is the case, in particular, when anisotropic 
flow is present, as an explicit calculation in Sec. Ill C will 
show. 

Therefore, the contribution of collective effects to the 
cumulants becomes dominant as k increases. Physically, 
the cumulant of order k isolates the genuine fc-particle 
correlation if k <C M: taking the logarithm in Eq. (15) 
effectively subtracts out the contributions of lower-order 
correlations, as shown explicitly in Rcfs. [33, 34]. In order 
to disentangle collective motion (which involves by defi- 
nition a large fraction of the particles) from lower-order 
correlations, it is therefore natural to construct cumu- 
lants of order k as large as possible. 

The idea of the cumulant expansion proposed in 
Rcfs. [33, 34] was to construct explicitly the cumulants 
at a given order k. Here, we propose to study directly 
the asymptotic limit when k goes to infinity, as we now 
explain. 11 The asymptotic behavior of the cumulants 
defined by Eq. (15), when k goes to infinity, is determined 
by the radius of convergence of the power series expan- 
sion, i.e., by the singularities of In G 6 {z) which are closest 
to the origin in the complex plane of the variable z. It is 
obvious from the definition, Eq. (6), that G (z) has no 
singularity. Hence the only singularities of h\G e (z) are 
the zeroes of G e (z). 

Let us denote by Zq the zero of G e (z) which is closest 
to the origin in the upper half of the complex plane. Re- 
member that Eq. (8) relates the behavior of G 6 (z) in the 
lower half of the complex plane to that in the upper half, 
so that we only need to study the latter. The asymp- 
totic behavior of the cumulants for large k is derived in 
Appendix C, Eq. (C3): 

a ««m> ~Hw)- (16 » 

As expected, large-order cumulants depend only on Zq. 



These subevents, which may contain only a few particles, have no 
relation with the "subevents" used to determine the event-plane 
resolution in the standard flow-analysis method. 



11 When k becomes as large as the multiplicity M, the cumulant no 
longer corresponds to the genuine fc-particle correlation. While 
the latter cannot be defined for k larger than M, the cumulants 
arc well defined to all orders. 
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Collective effects, if any, are uniquely determined by 
Zq. They result in larger correlations, i.e., in higher val- 
ues of the cumulants than in the absence of collective mo- 
tion. According to Eq. (16), this means that Zq should 
be smaller if collective effects are present. In particular, 
z^ will in general come closer to the origin as the size of 
the system, M, increases, as we shall see in Sec. IIIC. If 
only short-range correlations are present in the system, 
on the other hand, Zq does not depend on M in the limit 
of large M. This will be shown by means of an explicit 
example in Sec. IV. 



C. Relation with anisotropic flow 

We now evaluate the generating function G s (z) in the 
presence of anisotropic flow. We assume that \z\ is much 
smaller than unity (with weights of order unity), which 
will be justified later in Sec. IV. Since all coefficients 
in the power-series expansion of Eq. (14) are of the same 
order of magnitude, we can truncate this series for \z\ <C 1 
by keeping only the first two terms, which we rewrite in 
the form 



\n(e z ^\^ R )^(Q e \^ R )z + 



a 2 z 2 



(17) 



In this equation, (Q 6 \^r) denotes the average value of 
Q e for a given Using the definition of integrated 

flow, Eq. (4), and symmetry with respect to the reac- 
tion plane, one obtains (Q x \$r) = V n cos(n$^) and 
(Qv\®r) = V„sin(n$ij). From the definition of Q 6 , 
Eq. (3), it follows that 



(Q 9 \^r) = V n cos(n(* fl - 0)). 



(18) 



The parameter a in Eq. (17) is the standard deviation of 
Q for a fixed <fr R : 



a 2 EE2(((Q e ) 2 |^)-(Q^) 2 ). 



(19) 



Comparing Eqs. (14) and (17), a 2 scales linearly with the 
size of the system. For independent particles, in partic- 
ular (no correlations, no anisotropic flow), using Eq. (3), 
one obtains 



a 




(20) 



From the average value of e z ® e for a fixed &r, given 
by Eq. (17), one deduces the probability distribution of 
Gf (for a fixed $^) by inverse Laplace transform. It 
is easily found to be Gaussian: this shows that the ap- 
proximation made in keeping only the first two terms in 
expansion (14) is the central limit approximation already 
used in Rcfs. [2, 27, 43]. 

With the help of Eq. (18), Eq. (17) can be rewritten 
in the form 



,zQ e 



<jj fl \ _ e <T 2 z 2 /4+V n cos(n(<S> a - 



8))z 



(21) 



We now average over <&r, so as to obtain G e (z). Here, 
we further assume that a in Eq. (19) is independent of 
§r. The validity of this assumption will be discussed 
in Sec. VA. The theoretical estimate of G e (z) under 
this condition is denoted by G c .\.(z), where the subscript 
"c.l." refers to the central limit approximation: 



G c± (z) = e a2z2/4 Io(V n z), 



(22) 



where Iq is the modified Bessel function of order 0. The 
generating function is independent of 9, as expected from 
azimuthal symmetry. It is even [a natural consequence 
of Eq. (7) and azimuthal symmetry], so that cumulants 
of odd order vanish. It is also an even function of V n , 
so that the sign of V n cannot be determined from the 
generating function. From now on, we assume that it is 
positive. 

Anisotropic flow enters the generating function 
through the combination V n z. Therefore, cumulants of 
even orders 2k are of order (V n ) 2k , i.e., they scale with 
M like M 2k as expected from the general discussion in 
Sec. IIIB. From the cumulant of order 2k, one can ob- 
tain an estimate of the flow V^{2k} in the following way: 
take the logarithm of Eq. (22), expand to order 2k, and 
identify the coefficient with the corresponding coefficient 
in Eq. (15). This is roughly the procedure proposed in 
Refs. [33, 34] (a more thorough comparison is performed 
in Appendix B). Here, we want instead to take directly 
the large-order limit k — > oo. The large-order expansion 
of In Gc.i. (z) is given by Eq. (C3), where zo denotes the 
"first" zero of G c .i. (z). All zeroes of Iq lie on the imagi- 
nary axis. On the imaginary axis, we have 



G c .Uir)=e- 2r2 ^J (V n r), 



(23) 



where Jo is the Bessel function of order zero, and r is 
real. The first zero of J lies at joi — 2.405, hence the 
first zero of G c .i.(z) is at 



z = ir 



mi 



(24) 



Since V n scales like the total multiplicity M, zq comes 
closer to the origin as M increases, as expected from the 
discussion in Sec. IIIB. Note that Gei.(z) has no zero 
when there is no anisotropic flow. 

Identifying the terms in z 2k in the expansions of 
In Gc.i. (z) and In G (z), one obtains for large k 



V"{2kY k = (-l) k { ]m ) 2k Ke 



1 



8\2k 



(4) 



(25) 



If the zero of G e (z) lies exactly on the imaginary axis at 
Zq = ir®, as the zero of G c .i.(z), then V®{2k} converges 
to a limit for large k, which is V„ {oo} defined by Eq. (9). 

If, however, Zq has a (small) non-vanishing real part, 
V„ {2k} does not converge for large k. Unfortunately, this 
is the general case: experimentally, the available number 
of events iV ev t s is finite, and the resulting statistical fluc- 
tuations of G e (z) push the zeroes of G 6 (z) slightly off 
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the imaginary axis. These deviations, however, are phys- 
ically irrelevant. This is why we suggested, in Sec. II B, 
to find the first minimum of |G e (zr)| (denoted by Tq), 
rather than the first zero of G®{z) in the complex plane. 
Both procedures are of course equivalent when z® lies on 
the imaginary axis, as it should in an ideal experiment. 
We then approximate Zq ~ ir® in Eq. (25). Then, we ob- 
tain a consistent limit for large k, which is again Eq. (9). 
However, one should keep in mind that V®{oo} is the 
limit of V®{2k} for large k only when the first zero lies 
exactly on the imaginary axis. 

The careful reader will have noted that the theoret- 
ical estimate G c .\.(ir) is real [see Eq. (23)]. One could 
then argue that the imaginary part of G e (ir) is also ir- 
relevant, and look for the first zero of the real part of 
G (ir), rather than the first minimum of |G 8 (ir)|. How- 
ever, this is incorrect. We shall see in Sec. VIII that 
small inhomogeneities in the detector acceptance slightly 
shift the phase of G e (ir). The real part oscillates and 
has zeroes, but they are irrelevant. 

The following Sections of this paper are essentially dis- 
cussions of the errors which occur when one of the hy- 
potheses listed in Sec. Ill A is violated. First, we dis- 
cuss errors resulting from the simplifications made in 
Sec. IIIC. They will be shown to amount to finite mul- 
tiplicity corrections [violations of hypothesis A)]. These 
errors are of two types: 1. Even if no anisotropic flow is 
present in the system, the generating function does have 
zeroes, so that the analysis will yield a spurious flow. Its 
magnitude is discussed in Sec. IV. 2. When anisotropic 
flow is present, there are finite multiplicity corrections to 
the formula (9). These are discussed in Sec. V. Next, we 
discuss a specific violation of hypothesis B), namely fluc- 
tuations of impact parameter in the sample of events, in 
Sec. VI. Violations of hypothesis C) result in statistical 
errors, which are computed in Sec. VII. Finally, viola- 
tions of hypothesis D), i.e., detector effects, are discussed 
in Sec. VIII. 

D. Relation with the Lee- Yang theory of phase 
transitions 

Lee and Yang showed in 1952 that phase transitions 
can be characterized by the distribution of the zeroes of 
the grand partition function in the complex plane [39]. 
The grand partition function is defined as 

= £ Z N e^ kT , (26) 
jv=o 

where Zj^ is the canonical partition function for N par- 
ticles at temperature T in a finite volume V. Both T 
and V are fixed. In order to make the analogy with our 
formalism more transparent, we rewrite the above grand 
partition function in the form 



where z = (/j, - H c )/kT, P N = Z N e^ N ' kT /Q{n c ), and 
/i c is a reference value. Lee and Yang rewrote the grand 
partition function in a similar way, but they choose the 
variable y — e z (fugacity) instead of z. 

Physically, the coefficients Pjy in Eq. (27) represent 
the probability of having N particles in the system at 
chemical potential \i c . The grand partition function can 
therefore be rewritten as a statistical average with this 
probability distribution: 

G{z) = (e zN ). (28) 

The formal analogy with our generating function, Eq. (6), 
is obvious. 

Lee and Yang studied the repartition of the zeroes of 
G(z) in the plane of the complex variable z. 12 There 
is no zero on the real axis, since Eq. (27) is a sum of 
positive terms. Lee and Yang first showed that if a phase 
transition occurs at ji = fi c , the zeroes of G(z) in the 
complex plane of the variable z come closer and closer 
to the origin z — as the volume of the system, V, 
increases [39]. If no phase transition occurs, on the other 
hand, zeroes remain at a finite distance from the origin. 

This approach can easily be extended to the canonical 
partition function, when expressed as a function of the 
variable z = l/(kT c ) — l/(kT), where T c is the critical 
temperature. In this case, one usually speaks of Fisher 
zeroes [49] rather than Lee- Yang zeroes. 

The properties of G (z) derived in Sec. IIIB also ap- 
ply to G(z) as defined by Eq. (28), provided one replaces 
Q e by the number of particles, N, and the multiplicity 
M by the volume V. More specifically, it is well known 
that if a system can be decomposed into n independent 
subsystems, the partition function can be factorized into 
the product of the individual partition functions of each 
subsystem: G(z) — YYj=i^j( z )- Then, the zeroes of 
G(z) for the whole system are the zeroes of the partition 
functions Gj(z) for each subsystem. If the subsystems 
are equivalent, the position of the zeroes of G(z) does 
not change as the number of subsystems increases. In 
particular, their distance from the origin does not de- 
crease as the system size increases. This property can be 
generalized to systems with short-range correlations. 

If long-range correlations are present, on the other 
hand, large-order cumulants are larger (see the discus- 
sion in Sec. IIIB), and the first zero of G(z) comes closer 
and closer to the origin as the system size increases. This 
can be easily understood in the case of a first-order phase 
transition, say, a liquid-gas phase transition. At the crit- 
ical point /j, — He, for a large system, one can have any 



These zeroes completely characterize the grand partition func- 
tion if there is a hard-core interaction: in this case, there is an 
upper bound on the multiplicity N for a finite volume, so that 
G(z) is a polynomial of the variable y = e z , which is completely 
determined by its roots and its value at the origin G(0) = 1. 
However, most of the crucial results of Lee and Yang are still 
valid if this hard-core assumption is relaxed. 
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mixture of the (low-density) gas phase and the (high- 
density) liquid phase. The probability distribution Pjy 
in Eq. (27), instead of being sharply peaked around its 
average value, is widely spread between two values iV m ; n 
(pure gas) and N max (pure liquid) which both scale like 
the volume V. Then, the partition function G(z) depends 
on the volume V essentially through the combination zV, 
and consequently its zeroes scale with the volume like 
1/V. Anisotropic flow produces a similar phenomenon: 
its contribution to the generating function G (z) is a mul- 
tiplicative term I (V n z) in Eq. (22), where V n scales like 
the total multiplicity M, and zeroes accordingly scale like 
1/M. Anisotropic flow appears formally equivalent to a 
first-order phase transition in this respect. The impor- 
tant difference with statistical physics is that the system 
size is much smaller in practice. As a consequence, ze- 
roes never come very close to the origin, but the physics 
involved is essentially the same. 

The analogy with Lee- Yang theory goes one step fur- 
ther. In a second paper [40], Lee and Yang showed that 
for a very general class of interactions, the zeroes of G(z) 
are located on the imaginary axis (i.e., on the unit circle 
for their variable y = e z ). Here, although we do not have 
a general proof for this result, the same property holds in 
most cases. In particular, we have seen in Sec. Ill C that 
zeroes resulting from anisotropic flow lie on the imagi- 
nary axis [see Eq. (23)]. The main reason is that our 
generating function G e (z) is even (for a perfect detector 
and an infinite number of events). Now, we know that 
zeroes come in conjugate pairs due to Eq. (8), i.e., if z is 
a zero, then Zq is also a zero. If z$ lies on the imaginary 
axis, Zq = — zq, which satisfies the parity requirement. 
This parity argument does not show by itself that zeroes 
should be on the imaginary axis, but it is nevertheless 
crucial in the proof by Lee and Yang, and it is interest- 
ing to note that their result remains valid in the cases of 
interest in the present paper. 

Finally, we wish to mention that Lee- Yang zeroes 
were also used in high-energy physics in another con- 
text, namely in analyzing multiplicity distributions. It 
was found that the generating function of multiplicity 
distributions has zeroes that tend to lie on the unit circle 
in the complex plane of the fugacity, as Lee- Yang ze- 
roes [50]. However, it was later shown that this behavior 
merely reflects general, well-known features of the mul- 
tiplicity distribution [51], so that this approach does not 
seem to give new insight on the reaction dynamics. 



E. Differential flow 

Let us now justify the recipes given in Sec. II C to an- 
alyze differential flow. We shall proceed in the same way 
as for integrated flow. We start with a discussion of our 
hypotheses and of their implications. Next, we define 
the cumulants of the correlations between an individual 
particle, whose differential flow we are interested in, and 
the flow vector Q e . We discuss the order of magnitude of 



these cumulants in the case when only short-range corre- 
lations are present in the system, and derive the general 
expression of their large-order behavior. Then, we com- 
pute their value in the presence of anisotropic flow, in 
order to relate them to the differential flow. Finally, we 
explain why "autocorrelations" arc negligible in our ap- 
proach. 



1. Preliminaries 

Our hypotheses are the same as in Sec. Ill A. Here, we 
want to analyze the flow in a given phase space window. 
We denote by the azimuthal of a particle belonging to 
the window (which we call a proton for sake of brevity) . 
In order to study its flow, we have to correlate it to the 
flow vector Q e , or equivalently, to the generating func- 
tion e zQ . In order to study the Fourier harmonic v' mn , it 
is natural to construct averages over all protons such as 
(cos(mn(V> — 0)) e z ® ). We can first perform this average 
for a fixed orientation of the reaction plane, Hypoth- 
esis B) allows us to obtain an equation similar to Eq. (13): 
the contributions of independent subevents factorize. Di- 
viding by (e z ® |$r), these contributions cancel, except 
for the subevent containing ip. One thus obtains an equa- 
tion similar to Eq. (14) 



cos(mn(?/> — 6)) e z< ^" |$k 



a' + b'z + c'z 2 • • • , (29) 



where the coefficients a', b' , d are independent of the 
system size (since they only involve the subevent to which 
i\) belongs), and typically of order unity [if weights in 
Eq. (3) arc of order unity]. This is the mathematical 
formulation of hypothesis B). 



2. Cumulants 

We first introduce the following generating function 
D e m (z) = {cos(mn(?/> - 6)) e zQ " } . (30) 
This generating function has the symmetry properties 
D 9 +*' n {z) = {-l) m Dl{-z) (31) 

and 

[Dl(z)\* = Dl{z*), (32) 

which are analogous to Eqs. (7) and (8), respectively 
Furthermore, if the number of events is infinite and if 
the detector has perfect azimuthal symmetry, D e m {z) is 
independent of 6. 

The cumulants arc defined by 



D e m {z) 



GHz) k\ 

v ; fc=0 



(33) 
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where G e (z) is defined by Eq. (6). If there is no 
anisotropic flow, outgoing particles are not correlated 
with the orientation of the reaction plane, and the lhs 
of Eq. (29) is the lhs of Eq. (33). This shows that cumu- 
lants are independent of the system size (and typically of 
order unity) when there is no anisotropic flow. 

When collective motion is present, on the other hand, 
cumulants scale like M k , as we shall see below in the 
case of anisotropic flow. Following the same reasoning as 
for integrated flow, this shows that the best sensitivity 
to collective flow is achieved by constructing large-order 
cumulants. 

The large-order behavior of the cumulants is deter- 
mined by the singularities of the generating function in 
the lhs of Eq. (33). Since D e m {z) is analytic in the whole 
complex plane, the singularities are here again the zeroes 
of G e (z). The approximate expression of the cumulants 
is derived in Appendix C. It is given by Eq. (C4): 



h dk 



Re - 



2 KM 

M fc+1 {Q e e*°Q e } 



(34) 



where zq is the first zero of G (z), and where we have 
evaluated the derivative of G e {z) in the denominator, 
using definition (6). 



3. Relation with anisotropic flow 

Let us now compute Df n (z) when there is anisotropic 
flow. We use Eq. (29) and keep only the first term in the 
power series expansion: 



cos(mn(ip — #)) e z ® \®i 



(cos(mn(V> - 8))\§r) . 



This amounts to assuming that the proton is not corre- 
lated with the other particles for a fixed orientation of 
the reaction plane, and that there is no autocorrelation, 
i.e., that the proton is given zero weight in the definition 
of the flow vector, Eq. (2). 

The denominator of the lhs of Eq. (35) is given by 
Eq. (21), while the rhs is given by a relation analogous 
to Eq. (18): 

(cos(mn(^ - 9))\$r) = v' mn cos(mn($ fl - 6)). (36) 

We denote by D mc ,\\z) the value of D e m {z) under the 
above hypotheses. Using the last two equations and av- 
eraging over $_r, one obtains 



D mc . l .^)^e a2z2 ^I m (V n z)v' m 



(37) 



Combining this identity with Eq. (22), we finally obtain 

(38) 



D mc .\.{z) I m {V n z) , 



The cumulants are obtained by expanding this function 
in powers of z, as in Eq. (33). The only non- vanishing 
terms are those of order z 2k+m , where k is a positive 
integer. 

One can thus obtain an estimate of differential flow 
from the cumulant of order 2k + m by expanding both 
Eqs. (33) and (38) and identifying the coefficients of 
z 2k+m j n t nese expansions. This requires to know the 
integrated flow V n , which has been estimated earlier. 

The large-order coefficients in the expansion of Eq. (38) 
are given by Eq. (C4), where D e m {z) and G (z) are re- 
placed by I m (V n z) and Io(V n z), respectively. The deriva- 
tive of I (V n z) with respect to z is V n I\(V n z), so that 



<h R / -2 I m {V n zo) 

k\ ~ G Uzo) fc+1 V n h(V n Zo) 



(39) 



Next, we replace V n by its estimate Eq. (9). The pole 
then lies at zq — ir®. Inserting this value into the previ- 
ous equation, and using the relation I m (ir) — i m J m (r), 
we obtain 



dk ( -2i m ~ 1 

k\ VK) fe+1 



Jmjjoi 

JiUoi)J v^y 



(40) 



This is nonvanishing only if k — Ik! + to, as expected. 

In order to obtain an estimate of v' mn , we must iden- 
tify this expression with the large-order expansion of the 
measured cumulants, given by Eq. (34). As explained 
above in the case of integrated flow, the real part of the 
first zero, z , is irrelevant, so we replace z in Eq. (34) 
by irfj. One thus immediately recovers Eq. (12). The 
remarkable result is that the generating function for dif- 
ferential flow D 6 m {z) need only be evaluated for a single 
value of z. In this respect, the analysis is simpler than 
the cumulant method of Refs. [33, 34], and more reliable 
numerically. 



4- Autocorrelations 

In the standard event-plane method [26], differential 
flow is obtained by correlating the particle of interest 
with the flow vector of the event, Eq. (2). It is necessary 
to first remove the particle from the definition of the 
flow vector, otherwise the resulting autocorrelation would 
produce a spurious differential flow. 

One must be aware that even after autocorrelations 
have been removed, an error of the same order of mag- 
nitude may still be present due to nonflow correlations: 
consider the simple example where particles are emitted 
in collinear pairs. 13 When one correlates a particle with 



G cl .{z) I (V n z) 



13 Such collinear particles are expected from minijets, as will be 
seen in Sec. IV. Decay products from high transverse momentum 
resonances will also be almost collinear. 
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the flow vector, if the flow vector involves the other parti- 
cle in the pair, the resulting correlation will be exactly of 
the same magnitude as an autocorrelation. Conversely, 
methods which are less biased by nonflow correlations, 
such as higher-order cumulants, are also less biased by 
autocorrelations [33]. 

Here, autocorrelations do not produce any spurious 
flow by themselves. To show this, we consider for sim- 
plicity a particle which is not correlated with any other 
particle, whatever the reaction-plane orientation, so that 
its differential flow vanishes. We are going to check that 
the above procedure yields the correct result v' mn = 0. 

Let us denote by Q' e the value of Q e [Eq. (3)] after 
subtracting the contribution of the particle under study: 
Q'° = Q e — wcos(n(ip — 9)), where w is the weight as- 
sociated with the particle. If the particle with angle ip 
is not correlated with the other particles, the average in 
Eq. (30) can be factorized: 



D e m {z) = (cos(mn(ip - 9)) e 



zw cos(n(ip— 0)) \ I zQ 1 ' 



(41) 

The last term in the rhs can be estimated as in Sec. Ill C. 
This gives an equation similar to Eq. (22): 



e*Q' e )=e°' 2 > 2 /*I (V n z). 



(42) 



The crucial point is that the particle which has been 
removed does not flow by hypothesis, so that it does 
not contribute to the integrated flow in Eq. (4). There- 
fore, Q e and Q' e yield the same integrated flow value 
V n , although their statistical fluctuations differ in gen- 
eral (a' < cr). 



Therefore. 



vanishes for the same value of z as 



G 6 {z). Now, our estimate of differential flow, Eq. (12), 
involves the generating function D e m {z) precisely at the 
point zq where G e (z) vanishes. This shows that no spu- 
rious differential flow appears, although autocorrelations 
have not been removed. 

On the other hand, when there is differential flow, au- 
tocorrelations produce spurious higher harmonics of the 
flow. This will be discussed in Sec. VB. 



IV. SENSITIVITY OF THE METHOD 

Even if there is no flow in the system, the method pre- 
sented in this paper will give a spurious "flow" value. In 
this Section, we estimate the order of magnitude of this 
spurious value, and show that it is smaller than with any 
other method of flow analysis. We assume that hypothe- 
ses C) and D) are satisfied, i.e., we neglect statistical 
fluctuations and detector effects. 

Since no spurious flow appears in the central limit ap- 
proximation, as shown in Sec. Ill C, we need to go beyond 
this approximation. In this Section, we choose to work 
out exactly a simple model, but the order of magnitude 
of our results is general. 



The model is the following: we assume that all par- 
ticles are emitted in jets, each jet containing q collincar 
particles, where q is not much larger than unity. Since 
the total multiplicity is M, the total number of jets is 
M/q. The case q = 1 corresponds to independent parti- 
cle emission. With unit weights, Eq. (3) becomes 



M/q 

Q e = q^cos(n( ( p j -e)), 



(43) 



where the 4>j are the azimuthal angles of the jets. 

We assume that the jets are mutually independent and 
uniformly distributed in azimuth. Equation (6) then 
gives 



G°(z) = [I (qz)] 



M/q 



(44) 



The generating function is even and independent of 9, 
thanks to the azimuthal symmetry. Note that In G 9 {z) is 
proportional to M, as expected from the general discus- 
sion in Sec. Ill B. 

As explained in Sec. Ill C, the estimate of V n from the 
cumulant of order 2k, denoted by V®{2k}, is obtained by 
expanding In G e (z) to order z 2k and identifying with the 
expansion of In G c . \.(z) [Eq. (22)] to the same order. One 
thus obtains 



V?{2fc} = q 



l/2fc 



(45) 



This spurious flow increases with q, i.e., with the mag- 
nitude of nonflow correlations. There is even a spurious 
flow for q = 1, i.e., for independent particles. This is 
a consequence of "autocorrelations" which appear when 
expanding G e (z) in powers of z, and are not present in 
other methods of flow analysis. An alternative choice for 
the generating function, which is free from these autocor- 
relations, is discussed in Appendix A. It will be shown 
that it does not give better results as soon as q > 2, i.e., 
as soon as nonflow correlations are present: when this is 
the case, Eq. (45) still gives the order of magnitude of 
the spurious flow. 

As expected from the general arguments in Sec. IIIB, 
the spurious flow decreases as the order k of the cumulant 
expansion increases. We now evaluate V% {oo} defined by 
Eq. (9). Although there is no flow in the system, the gen- 
erating function (44) does have zeroes on the imaginary 
axis. The first one lies at z = ir^, with t-q = jo\/q. [Note 
that it is a multiple zero, because of the power M/q, so 
that the rhs of Eq. (16) must be multiplied by M/q] In 
agreement with the general discussion in Sec. IIIB, the 
position of the zero does not depend on the multiplicity 
M. Equation (9) then gives a spurious "flow" value 



V£{00} = q. 



(46) 



It coincides with the limit of V£{2k} [Eq. (45)] for large 
k, as expected. The order of magnitude of this result is 
general if weights are of order unity. 
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Flow can be identified unambiguously only if it much 
larger than the spurious value. Using Eq. (5), and the 
assumption that q is not much larger than unity, one 
finally obtains the condition under which our method can 
be applied: 



Vn>> M- 



(47) 



By comparison, the conditions under which one can ap- 
ply the standard flow analysis (or two-particle meth- 
ods) or four-particle cumulants are v n 3> 1/M 1 / 2 and 
v n 3> 1/M 3 / 4 , respectively. The present method is more 
sensitive, in the sense that its validity extends down to 
smaller values of the flow. 

The condition (47) implies that the first zero of G e (z), 
approximately given by Eq. (24), satisfies the condition 
\zo | <C 1, which justifies the approximation made at the 
beginning of Sec. Ill C. As already discussed in Refs. [33, 
34] , it is probably impossible to measure anisotropic flow 
if this condition is not satisfied. 

Note that the present model of nonflow correlations 
is a very extreme one, since we have assumed that they 
affect all particles, and that particles within a jet are ex- 
actly collincar. More realistic estimates of the spurious 
flow would give significantly smaller values. However, re- 
fining the above estimates would be academic: we shall 
see in Sec. VII D that if there is no anisotropic flow in 
the system, the spurious flow created by statistical fluc- 
tuations is much larger than the flow created by nonflow 
correlations, unless the number of events is unrealistically 
large, so that the present method will not be limited by 
nonflow correlations in practice. 



V. SYSTEMATIC ERRORS 

In this Section, we estimate the order of magnitude 
of the systematic error due to the method itself on our 
flow estimate. We discuss first integrated flow, then dif- 
ferential flow. We assume that hypotheses C) and D) 
are satisfied, i.e., we neglect statistical fluctuations and 
detector effects. 



A. Integrated flow 

The error on the integrated flow is a consequence of 
the two approximations made in Sec. Ill C. The first one 
was to keep only the first two terms in the power-series 
expansion of Eq. (14). The second one was our assuming 
that a in Eq. (19) was independent of 

We first estimate the error resulting from the term 
proportional to z 3 in Eq. (14), which was neglected 
in Eq. (17). We know from the general discussion in 
Sec. Ill A that the coefficient c in front of z 3 is indepen- 
dent of the size of the system, M, and vanishes when 
there is no anisotropic flow. Hence, it is natural to as- 
sume that it is of order v n . Since the value of z of interest 



is of order l/(Mv n ) [from Eq. (24)], the correction to the 
rhs of Eq. (17) arising from the z 3 term is of the order 
of Mcz 3 ~ l/(Mv n ) 2 . Since the lhs of Eq. (17) is a 
logarithm, this is a relative correction to the generating 
function. Hence, the relative error on the flow estimate 
will be of the same order of magnitude. 

Next, we estimate the error due to the dependence of 
(T on <f>/f. A detailed calculation of this effect is per- 
formed in Appendix D. In particular, it is shown that 
anisotropic flow results in contributions to a 2 propor- 
tional to sin 2 (n(<E>fl — 9)). These contributions consist 
of two terms, which are of magnitude Mv 2 and Mv2n 
(i.e., involving a higher harmonic), respectively. Now, 
a 2 is multiplied by z 2 in Eq. (17). Replacing z by the 
value of interest in Eq. (24), ^-dependent terms yield 
contributions to the rhs of Eq. (17) of orders 1/M and 
i)2n/(M^), respectively. As explained above, the rela- 
tive error on the flow estimate due to these terms will be 
of the same order of magnitude. 

Finally, note that the contribution of the term of order 
z 4 to the rhs of Eq. (14) may be of the same order of 
magnitude, or even larger, than the term of order z 3 , 
although |z| is much smaller than unity. This is because 
the coefficient in front of z 3 is of order v n , which can 
be much smaller than unity, while the coefficient in front 
of z 4 is of order unity. However, only the part of this 
coefficient which depends on matters: indeed, the 
<I> i^-independent part gives a multiplicative contribution 
to G (z) of the type exp(Az 4 ), which has no zeroes. Now, 
the ^-dependent contribution to the term of order z 4 
is much smaller than the <& /{-contribution to the term of 
order z 2 , which has been estimated above. 

In summary, the order of magnitude of the relative 
systematic error on the flow is 



V n {00} - V n 



+ O 



o 



1 



{Mv n f 



V2n 

Mv 2 



(48) 



All three terms involve powers of 1/M: they can be 
viewed as finite multiplicity corrections, resulting from 
the violation of hypothesis A) in Sec. Ill A. 

Let us comment on the magnitude of these errors. The 
first term, a relative correction of order 1/M, is always 
small since M ^> 1. The last two terms are much less 
trivial. They also appear when flow is analyzed from the 
cumulant of four-particle correlations, where they result 
from an interference between correlations due to flow and 
nonflow correlations, as discussed in detail in Appendix 
A in Ref. [33]. As long as condition (47) is satisfied, the 
term of order l/(Mv n ) 2 is a small correction. The rel- 
ative magnitude of the first two error terms in Eq. (48) 
actually depends on the value of the resolution parame- 
ter x which will be defined below in Sec. VII A, and is 
of order v n \fM. The largest correction term in Eq. (48) 
is the first one for large x, and the second one for small 
X- Higher harmonics are generally of smaller magnitude: 
typically, V2n ~ v n , so that the third term in the rhs of 
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Eq. (48) is comparable to the first term. An important 
exception is that directed flow is much smaller than ellip- 
tic flow at ultrarelativistic energies. In this regime, the 
third term dominates the error for n = 1. 

Nonflow correlations involving a few particles increase 
the systematic error on the flow, but the order of magni- 
tude is still given by Eq. (48). In the extreme case where 
all particles are emitted in bunches of q collinear parti- 
cles, as in Sec. IV, one must replace in Eq. (48) M by 
M/q, which is the effective multiplicity, i.e., the number 
of independent angles (the number of independent de- 
grees of freedom). Therefore the first and third term in 
Eq. (48) will be increased by a factor of q, and the second 
term by a factor of q 2 . 

Finally, how does our systematic error compare with 
the systematic error on the flow estimate V n {2k} from 
cumulants of finite order 2k? Possible nonflow correla- 
tions between 2k particles yield an additional term: 14 



V n {2k} - V n 



' M 



O 



O 



1 



{Mv n f 

V2 ^\ + o 



Mvl) 



1 



J^2k-l y 2k 



• (49) 



In particular, this equation can be used to show that 
the error on flow estimates from standard (event-plane 
or two-particle) methods is much larger than with the 
present method. Setting 2k = 2 in the previous equation, 
the last term dominates over the first three terms, so that 
we may write 



Vn{2} - V n 
V n 



= o 



1 



Mvl 



(50) 



With higher-order cumulants, 2k > 4, the last term in 
Eq. (49) dominates only if 



1 



Vn ^ n/r\-X/(2k-2') ' (51) 
tile error was given correctly in 



/ma" 



14 The magnitude of the sys* 

Appendix A of Rcf. [33] , but not in Ref. [34] where only the last 



If v n is larger than this value, the error is similar with the 
cumulant method and with the new method presented 
here. 



B. Differential flow 

The determination of differential flow relies on the pre- 
vious determination of integrated flow: one thus natu- 
rally expects a relative error on the differential flow of 
the same order as the relative error on the integrated 
flow, Eq. (48). This can be checked explicitly, see Ap- 
pendix D. 

An additional systematic error is due to the approxi- 
mation made in writing Eq. (35). It was assumed that 
the proton and the flow vector are independent. How- 
ever, they may be correlated, due to autocorrelations (if 
the proton is involved in the flow vector) or to nonflow 
correlations, whose effects are always similar to those of 
autocorrelations (see Sec. HIE 4). 

We now compute the error due to autocorrelations. For 
this purpose, we compute the next-to-leading term b'z in 
the power-series expansion Eq. (29). We compute only 
the contribution of the numerator of the lhs to this term. 
The contribution of the denominator yields a relative er- 
ror on v' mn of the same order as for the integrated flow, 
Eq. (48). If the proton has weight w' in the flow vector, 
a simple calculation gives 



b' 



w 
~2 



( v '( m -i)n cos[(m - l)n($R - 9)] 
+ v{ m+1)n cos[(m + l)n($ R -6)}). (52) 



Taking this term into account in Eq. (35) and the follow- 
ing equations, D 6 m {z) is no longer given by Eq. (37) but 

by 

term was kept, which was a mistake. 



J 



D° m (z) = e° 



I m {V n z)v' m 



2 7 " 



i{V n z)w'v[ m _ 1)n 
I 



-I m +i(V n z)w'v / 



(m+l)n 



(53) 



Since higher harmonics are generally of smaller magni- 
tude, we neglect the last term in the rhs. Evaluating the 
remaining two terms at z = iro = ijo\/V n , the differen- 
tial flow estimate v' mn {oo} is given by 



i{°°} 



+ 



joi./m-l(joi) W V (m-l)r. 



2Jm(joi) 



v„ 



(54) 



For m = 1, the extra correction term vanishes since 
■^o(ioi) = 0. One can show that the next term c'z 2 in 
the expansion Eq. (29) produces a relative error on v' n of 



the same order as that on the integrated flow, Eq. (48). 
For to = 2, on the other hand, we obtain 



'2n 



{oo} = v' 2n + 



JOI W ' V 'n 
4 V n ■ 



(55) 



Numerically with the input values M = 300 particles, 
v 2 — 6%, v' 2 — 7%, v' 4 = 3% used for the simulation of 
Sees. II B and II C, we find ^{oo} — 3.56%, in agreement 
with the result of the analysis when autocorrelations are 
not removed, <{oo} = 3.60 ± 0.29%, see Fig. 3. 



17 



As a conclusion, the systematic error on the differen- 
tial flow in the lowest harmonic v' n is a relative error of 
the same order of magnitude as for the integrated flow: 
it is small, and autocorrelations need not be removed. 
The situation is quite different for the higher harmonics 
v' 2n , v' 3n ..., where autocorrelations give a sizeable con- 
tribution. This contribution can be removed analytically 
using Eq. (55). However, one must be aware that an error 
of the same order of magnitude may remain as a result 
of nonflow correlations. The same conclusion holds with 
other methods of flow analysis. 

This error on higher harmonics has two consequences. 
At ultrarelativistic energies, where the flow vector is 
constructed with n = 2, the higher harmonic V4 is ex- 
pected to be of much smaller magnitude [52]. Since, 
furthermore, nonflow correlations are known to be size- 
able [18, 24], one may doubt whether V4 can be mea- 
sured at all. At lower energies, elliptic flow v 2 is usu- 
ally measured using the flow vector from directed flow 
n = 1 [5, 14, 44], this bias should be considered seri- 
ously. In particular, it may be large for mesons, which 
often originate from resonance decays and have therefore 
strong nonflow correlations. 



VI. FLOW FLUCTUATIONS 

As explained in Sec. Ill A, our results rely on the hy- 
pothesis that we are considering events with exactly the 
same impact parameter. In practice, however, the set of 
events in a given centrality bin spans some impact pa- 
rameter interval. In this Section, we are going to study 
the influence of such impact parameter fluctuations on 
the flow analysis [24]. 

We assume that the hypotheses of Sec. Ill A are sat- 
isfied for a fixed impact parameter b, and we denote the 
integrated flow at this impact parameter by V n {b). If b 
fluctuates in the sample of events, one can compute the 
generating function G e (ir) in two steps. One first av- 
erages over events with the same impact parameter b: 
(e* r< 3 \b) is given by the rhs of Eq. (23), where both V n 
and a may depend on b. Then, one averages over b: 

G e {ir) = (j (V n (b)r) e -(^ 2 /4^ , ( 56 ) 

where angular brackets denote an average over b. 

Our flow estimate V„{oo} is related to the first root, 
r 9 , of the equation G 9 {ir e ) = 0, by Eq. (9). First of 
all, if only a{b) fluctuates, while V n (b) is a constant, the 
position of the zero does not depend on b: this shows that 
the estimate of the flow is not affected by fluctuations of 
a if only those are present. Next, if V n (b) has small 
fluctuations around an average value, i.e., V n (b) = V n + 
SV n (b), then expanding Eq. (56) to first order in SV n (b), 
one easily shows that 14{oo} = V n , that is, the analysis 
using the present method yields the average value of the 
flow. 



When fluctuations are large, our estimate Ki{oo} docs 
not in general coincide with the average value of the flow, 
V n . It is sensitive to the fluctuations of V n (b), and also 
to the fluctuations of the width <r(b), which is involved in 
Eq. (56). However, it is possible to show under rather 
general conditions that if V n (b) lies between (V n ) m i n 
and (14) max , then 14{oo} also lies between (14) m ; n and 
(V n )max- For this purpose, we have to show that the 
first zero of G 9 (ir) defined by Eq. (56) occurs between 
t = joi/(K)max and r = joi/(V n ) min . Since J (x) > 
for < x < joi, and V n (b) < (V„) max , Eq. (56) shows 
that G e (ir) > for r < joi/(Ki)max- To prove our result, 
it is sufficient to show that G (ijoi/(Ki)min) < 0. Using 
Eq. (56), this holds as soon as Jo(joiV n (b) / (V n )mm) < 
for all b. Now, Jo (a;) is negative for x between j i and 
its second zero, j'02. Hence our result holds as soon as 



ioi < ioi tttt < .702 (57) 

i *n Jmin 

for all b. One obtains the condition: 



v .:;; max < ^ ~ 2.3. (58) 

(Vn)min JOI 

This condition is sufficient but not necessary. It shows 
that if the flow does not vary by a factor of more than 
2.3, the flow estimate V^,{co} lies between (V n ) m i n and 
(V n )maK, which is quite reasonable. 



VII. STATISTICAL ERRORS 

In this Section, we estimate the statistical errors which 
arise when the number of events N evts is finite. As with 
other methods of flow analysis, these statistical errors de- 
pend on the resolution parameter \- I n Sec. VII A, wc 
define this parameter, explain how it can be measured ex- 
perimentally, and give its approximate value for a variety 
of heavy-ion experiments. 

We then study, in Sec. VII C, the fluctuation of the 
generating function G e (ir), evaluated on the imaginary 
axis (i.e., for real r), around the true statistical average 
G c .i.(* r )- Then, in Sec. VII D, we show that even if there 
is no flow in the system, the data analysis through the 
present method will yield a non-zero "flow" value, which 
we estimate: we consider it as the lower bound on the 
flow, below which the present method cannot be used. 
When there is flow in the system, and when it is larger 
than this limiting value, statistical errors are small. They 
are estimated in Sec. VII E for the integrated flow and in 
Sec. VII F for the differential flow. 
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A. The resolution parameter 

An important parameter in the flow analysis is the 
resolution parameter \ defined as 15 



X 



Vn 

a 



(59) 



where V n and a are defined in Eqs. (4) and (19), respec- 
tively. Physically, it characterizes the relative strength of 
the flow compared to the finite-multiplicity fluctuations. 

With unit weights, V n = Mv n [Eq. (5)] and a = \[M 
for independent particles [Eq. (20)], so that \ — v n \[M. 
This shows that \ increases with both the flow and the 
number of detected particles. It is small for central colli- 
sions where v n is small due to azimuthal symmetry, and 
also for peripheral collisions where the multiplicity M is 
small. It is generally maximum for intermediate central- 
izes. Here again, we want to emphasize that one should 
use all detected particles in the analysis, so as to maxi- 
mize x- 

This parameter is related to the so-called "event- 
plane resolution" in the standard flow analysis. The 
event-plane resolution is defined as the average value of 
cos A$/f = cos(n(<!>„ — where n$> n is the azimuthal 

angle of the flow vector and $^ that of the reaction plane. 
It is related to \ by [27] : 



(cosA$ fl H^ xe -* 2 /2 



+ h 



(60) 



In Table I, we have summarized the values of x f° r 
various heavy-ion experiments in mid-central collisions 
(where x is largest), calculated with the help of Eq. (60) 
using the values of the event-plane resolution given in the 
cited references. These values of x ue typically between 
0.5 and 2.5. Note, however, that these values may not 
reflect only flow, but could be contaminated by nonflow 
effects. Such an instance, namely the influence of non- 
flow correlations from global momentum conservation on 
the determination of the first-harmonic event-plane res- 
olution, was discussed in Ref. [53]. 

How can x be determined experimentally? The numer- 
ator of Eq. (59) is the integrated flow V n , which is de- 
termined following the procedure introduced in Sec. II B. 
The denominator, a, is most simply obtained by averag- 
ing Eq. (19) over 9. In order to take into account possible 
azimuthal asymmetries in the detector acceptance, which 
were neglected in Sec. IIIC, we replace Eq. (18) with 

(Q 6 \^b) = V n cos(n(<5> R -9)) 

+ (Q x ) cos(n9) + (Q v ) sin(n0), (61) 

where we have used the first line of Eq. (3). The two 
additional terms vanish with a symmetric detector. One 



15 Please note that this definition of \ differs by a factor of l/\/2 
from the one adopted in Ref. [28] . 



TABLE I: Values of the resolution parameter \ an d the "dis- 
persion" factor used in the standard flow analysis for various 
experiments, n is the Fourier harmonic used to determine the 
event plane (1 for directed flow, 2 for elliptic flow). 



X 


(cos(A$ fl )) 


n 


System 


E/A 


Exp. 


Ref. 


0.90 


1/1.56 


1 


Au+Au 


90 MeV 


FOPI 


[6] 


2.6 


1/1.04 


1 


Au+Au 


400 MeV 


FOPI 


[6] 


1.65 


0.89 


1 


Au+Au 


2 GeV 


E895 


[12] 


0.52 


0.43 


1 


Au+Au 


8 GeV 


E895 


[12] 


0.39 


0.33 


1 


Pb+Pb 


158 GeV 


NA49 


[18] 


0.67 


0.53 


2 


Pb+Pb 


158 GeV 


NA49 


[18] 


1.25 


0.8 


2 


Au+Au 


65+65 GeV 


STAR 


[24] 


0.48 


0.4 


2 


Au+Au 


100+100 GeV 


PHENIX 


[21] 



then obtains 



(Ql + Ql)-(Q x ) 2 ~(Q v y 



(62) 



The value of a obtained in this way is more accurate than 
the approximate value in Eq. (20), which underestimates 
a when large nonflow correlations are present. Since sta- 
tistical errors are very sensitive to Xi as we shall see, it 
is important to estimate a as accurately as possible, in a 
model independent way, which is possible using Eq. (62). 



B. Preliminaries 

Here we derive simple results which will be useful in the 
remainder of the Section. The estimate of the integrated 
flow, V^{oo}, is obtained by averaging V^fjoo} over sev- 
eral values of 9, see Eq. (11), and a similar equation for 
differential flow: 



(63) 



fc=0 



In order to estimate the standard deviation of the aver- 
age P„{oo}, one needs to evaluate not only the standard 
deviation of P^{oo}, but also the correlation between dif- 
ferent estimates V®{oo} and V® {oo} with 6^6'. For 
this purpose, we shall have to evaluate averages such as 

(e ,r ^ e lr ® e ). We introduce the complex notation 



(64) 



One should note that the number z thus defined has noth- 
ing to do with the variable z involved in Eq. (6). With 
this notation, we have 



e irQ" _ e i(z*Q+zQ*)/2 



where Q = Q x + iQ y . Then, obviously, 

e irQ e e ir'Q e ' _ +z'*)Q+(z+z')Q*}/2 _ 



(65) 



(66) 
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We assume that the expectation value of e lrC2<> is given / e *rQ e e -ir'Q e \ — q c1 ^ z _ z >\^ (57) 

by our theoretical estimate, Eq. (23). Then, * ' 

rrQ e e ir'Q e '\ = Q cl .(i\z + z'\) from which we obtain 



Re(e lrQe )Re(e lr ' Qe )) = 



lm(e irQe )lm(e lr ' Q,> ) 



G cA .{i\z-z'\) + G cA {i\z + z'\) 
G cA .(i\z-z'\)-G cA .(i\z + z>\) 



(68) 



r 



All these identities will be used below. 



C. Fluctuations of G e (ir) 

In this Section, we estimate the magnitude of the 
statistical fluctuations of (the complex-valued) function 
G (ir) around its (real- valued) average G c .\.(ir). In 
particular, we derive the fluctuations of the modulus 
\G e (ir)\, which is the quantity minimized in Sec. II B 
to extract integrated flow. 

Let us first recall general results about sample aver- 
ages, i.e., average values evaluated from a finite sample 
of events. If x is an observable measured in an event 
(multiplicity, transverse energy, etc.), we denote by {x} 
the average value of x over the available sample of events, 
which is also called the sample average: 



Nevt 



{x} 



1 " ovts 



(69) 



The expectation value, which is the limit of this quantity 
when the number of events goes to infinity, is denoted by 
angular brackets as in the previous Sections. Note that 
(W) = (x). 

If iV cv t s is large but still finite, the sample average dif- 
fers from the expectation value by a small fluctuation 5x: 
{x} = (x) + Sx. If y denotes another observable associ- 
ated with each event, then the covariancc of the sample 
averages {x} and {y} is 

(Sxdy) = ({x}{y})-(x)(y) 

= ■^—((xy)-(x)(y)), (70) 

where in the second line we have used the property that 
events are statistically independent. Note that the sta- 
tistical properties of the fluctuations 5x, 5y are Gaussian 
due to the central limit theorem, hence they arc com- 
pletely determined by their covariancc: higher-order mo- 
ments are obtained using Wick's theorem. 

We can apply Eq. (70) to compute the statistical fluc- 
tuations of the generating function, Eq. (6), for z on the 



imaginary axis (z = ir, with r real). We introduce the 
decomposition 



G e (ir) = G c . l .(ir)+SG 9 (ir), 



(71) 



where SG (ir) denotes the fluctuation of G e (ir) around 
its expectation value. 

Let us first estimate the magnitude of 5G 6 (ir). For 
this purpose, we make use of Eq. (70) with x 
and y = e~ trQ : 



<I^V)I 2 ) 



1 



(72) 



At the zero of G c .i.(ir), the standard deviation of \G 6 (ir)\ 
is exactly l/^N cvts . This justifies Eq. (10). For large r, 
the expectation value G c .\.(ir) g° es to zero [see Eq. (23)], 
while the fluctuations do not: \SG e (ir)\ is of order 
l/VK^s- This can be easily understood: for large r, 
e lr ® is a random complex number uniformly distributed 
on the unit circle, so that the sample average G e (ir) in 
Eq. (3) is a random walk of iVevts steps, of length 1/N evts 
each, hence the order of magnitude of the fluctuations. 

We can now estimate the value of r for which the fluc- 
tuation 5G e (ir) becomes of the same order of magni- 
tude as the expectation value G c .i.(zr). The decrease of 
G c .\.(ir) for large r in Eq. (23) is dominated by the expo- 
nential factor, while Eq. (72) shows that \5G (ir)\ is of 
order 1 / y/N cv ts, so that both become comparable when 



1 



2 r 2 /4 



(73) 



cvts 



We denote by r c the corresponding value of r: 

r c = -^2\nN cvts . (74) 

G 

Fluctuations are relatively small below the critical value 
r c and dominate above. In the case of the Monte-Carlo 
simulation presented in Fig. 1, for instance, taking a = 
\[M, one obtains r c ~ 0.26: it is larger by a factor of 
2 than the position of the first minimum, Tq, which is 
the reason why the statistical fluctuations on our flow 
estimates, shown in Fig. 2, are small. 
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We finally estimate the fluctuations of the modulus 
| G s (ir) | , which is involved in determining the integrated 
flow (see Sec. II B). We assume that r < r c , so that 
|<5G e (?r)| is much smaller than G c .\.{ir). We separate 
the fluctuation into its real and imaginary part: 

G 9 (ir) = G c .i.(ir) +RcSG e (ir) +ilmSG e (ir). (75) 

Taking the modulus of Eq. (75), we then obtain, to lead- 
ing order in 8G e (ir) 



\G (ir)\ ~ iGd.(w-) + ReSG e (ir)\ 



i.e., only the real part of the fluctuations of G e (ir) con- 
tribute to the fluctuations of the modulus \G e (ir)\. Their 
covariance is obtained from Eq. (70), in which we take 
x = Re SG e (ir) and y — ReSG 6 {ir'), and with the help 
of Eq. (68): 



(76) 



Re SG 8 (ir) Re SG 9 ' {ir') 



2JV„. 



[G cA .(i\z - z'\) + G cA .(i\z + z'\) -2G cA .(ir) G cA .(ir')\ , 



(77) 



where we have used the property that the expectation 
value G cA .(ir) is real [see Eq. (23)]. 



D. Spurious flow from statistical fluctuations 

Let us now assume that there is no anisotropic flow 
in the system, V n = 0. In that case, analyzing the data 
with the present method will nevertheless yield a spuri- 
ous "flow" value, due to statistical fluctuations. In this 
Section, we first estimate the maximal such spurious flow 
value which analyses would give in the absence of real 
flow. This maximal value is in our eyes the minimal flow 
value which the method can safely reconstruct. Then we 
show that any V n extracted with the method larger than 
this value can in turn be reliably attributed to flow. 

If there is no flow in the system, Eq. (23) reads 



G cA .(ir) 



2 /4 



(78) 



so that there is no zero, nor any minimum of the true ex- 
pectation value G c .i. (ir) for finite r . With a finite number 
of events, however, |G e (-ir)| does in general have a min- 
imum, so that the procedure outlined in Sec. II B will 
yield a non-zero estimate of the flow, V®{oo}, which is 
unphysical. 

This minimum is due to the fluctuation SG e (ir): it 
typically occurs when SG e (ir) is of the same order of 
magnitude as the average value G c .i.(ir). The position 
of the minimum, r®, is therefore of order r® ~ r c , where 
r c is given by Eq. (74). The resulting "spurious flow" is 
given by Eq. (9). To this spurious flow, we can associate 
a resolution parameter as in Eq. (59) : 



A " 



(79) 



with V n = 0.] Using Eqs. (9) and (74), this yields 



Joi 



V2 In N c . 



(80) 



With iVe Vts = 20000 events, one obtains numerically 
X° = 0.54. If the analysis yields this (or a smaller) value 
for Xi the result cannot be attributed to flow, and is just 
an artifact due to statistical fluctuations. Experiments 
for which the resolution parameter \ ( see the values in 
Table I) is not larger than this spurious value will not 
have enough statistics to implement the present method, 
unless a proper choice of weights can significantly in- 
crease the value of \- 

The values of \ f° r actual heavy-ion experiments, listed 
in Table I, are often larger than this value, but not 
much larger, which might look somewhat worrying at 
first sight. Fortunately, as soon as the analysis yields a 
value slightly higher than x° m Eq. (80), the result can 
safely be attributed to collective flow, as we shall now 
show. 

As explained in Sec. VII C, a minimum of |G e («r)| oc- 
curs approximately when the real part of G (ir) van- 
ishes, i.e., when — ReSG e (ir) = G cA .(ir), sec Eq. (76). 
Now, ReSG e (ir) has Gaussian fluctuations, whose width 
is given by Eq. (77), in which we set z = z', and neglect 
the last two terms in the brackets in the rhs (which are 
at most of order l/y/N cv ts), so that 



Re5G e (ir)Y 



1 



2iVn 



(81) 



We can say with 98% confidence level that — ReSG e (ir) 
is smaller than twice this standard deviation, i.e., 



-Rc 5G e (ir) < yj 



(82) 



[Experimentally, a should be evaluated using Eq. (62) Therefore, if a minimum of |G (ir)\ occurs at r = 7q, we 
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can say with the same confidence level that 



GcJ.(*rg) < 



(83) 



Using Eqs. (9), (78) and (79), we obtain the following 
inequality with 98% confidence level: 



X" < 



Joi 



y/2 \n(N cvts /2) 



(84) 



With the same value N ev t s — 20000 as above, the upper 
bound is 0.56, only very slightly larger than the estimate 
from Eq. (80). This shows that even a value of the reso- 
lution parameter slightly larger than the rhs of Eq. (80) 
cannot be attributed to statistical fluctuations alone. In 
addition, this discussion holds for a fixed value of the 
reference angle 8. Averaging over several values of 8 de- 
creases the statistical error, and also slightly decreases 
the value of the spurious flow. This is why we gave 
X = 0.5 as the limiting value below which the method 
cannot be applied reliably in Sec. HE. 

In order to illustrate this discussion, we have simu- 
lated a data set with N cvts = 20000 events of multiplicity 
M = 300 each. Events were simulated with zero elliptic 
flow, V2 — for all particles, and we assumed a perfectly 
isotropic detector. The data set was then treated ac- 
cording to the procedure presented in Sec. II, using unit 
weights Wj — 1 and n = 2 in Eq. (2). The resulting 
V 2 e {oo}/M is displayed in Fig. 5 as a function of 8, to- 
gether with the upper bound given by Eqs. (79) and (84), 
in which we set a = \f~M. One first notes that all esti- 
mates V 2 e {oo} are smaller than the bound, which means 
at once that they should not be attributed to flow. The 
average over 8, Va{oo}/-M = 2.56%, is well below the 
upper bound. In addition, the values of V 2 {oo} strongly 
depend on 8, while we recall that in Fig. 2, also obtained 
in the case of a perfect detector, they were roughly inde- 
pendent of 8. 



Finally, a close look at the figure shows that the varia- 
tion of V 2 ' 1 {00} with 8 has four discontinuities, which can 
be qualitatively understood in the following way: when 8 
varies, G e (ir) varies continuously, but the first minimum 
of |G e (ir)| may disappear for some value of 8. Then, the 
second minimum becomes the first, and V 2 e {oo} suddenly 
jumps to a lower value (and vice- versa when a minimum 
appears below the first minimum). Such discontinuities 
are not expected when the minimum is due to anisotropic 
flow, and are by themselves a hint that the observed flow 
estimate is an effect of statistical fluctuations. 16 




FIG. 5: Analysis of N cvtB = 20000 events of multiplicity M = 
300, simulated without anisotropic flow. Crosses show the 
value of the spurious flow given by our analysis, V2{oo}/M, 
as a function of 9. The solid line displays the 98% CL upper 
bound on this spurious flow (see text). 



E. Statistical error on the integrated flow 

In this Section, we again assume that there is flow in 
the system, and that it is larger than the limit determined 
in the previous Section ("spurious flow"), below which 
the method cannot give reliable results. In that regime, 
we compute the statistical error on the integrated flow 
estimate V®{oo}, and on its average over directions 8, 
due to the limited number of events. We shall in particu- 
lar see that averaging several (typically 4 or 5) values of 
V„ {00} results in a decrease in the statistical uncertainty 
of a factor of almost 2. 

As explained in Sec. II B, integrated flow is determined 
from the first minimum of |G e (ir)|, which is denoted by 
7q. According to the discussion in Sec. VII C, for small 
fluctuations, this minimum is to leading order the first 
zero of Gc.i.(ir) + KeSG e (ir). It differs from ro [the first 
zero of G c .i.(ir), given by Eq. (24)] by 



ReSG e {ir) 



dG c .\.{ir) 




dr 


r=r 



(85) 



Using Eq. (23), and replacing ro with its value joi/V n 
[Eq. (24)], one obtains the denominator of this expres- 
sion: 



dG c .\.(ir) 
dr 



-V n Ji(joi) exp 



ioi 

4 X 2 



(86) 



The appearance of such discontinuities caused by statistical fluc- 
tuations, which do not occur for collisions with real flow, may 
be a good way to discriminate between both effects in the case 



of weak flow, slightly above the limit Eq. (80), provided one 
computes estimates V®{oo} for a large number of values of 6, 
instead of the 4 or 5 values which are enough to guarantee "low" 
statistical errors as argued in Sec. VII E. 
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where x is the resolution parameter denned by Eq. (59). 
Using now Eq. (9), and reinserting Eq. (86) in Eq. (85), 
one deduces the statistical fluctuation of the integrated 
flow estimate, 8V® = V® {oo} - V n : 

V n r 

= _ exp ( J oi/( 4 ^ 2 )) BjeSG e (ir ). (87) 
JoiJiUoi) 

The correlation (5V%5V„ ) can then be evaluated using 



Eq. (77). In this equation, the last term in the rhs van- 
ishes since r = r' = r , and G c .i.(ir ) = by definition 
of ro. In order to evaluate the first and second terms in 
the rhs, we use the expression of the expectation value 
G c .i.(ir), Eq. (23), and the two straightforward identities 
[see Eq. (64)] 

\z + z'\ = 2r cos(n(6-9')/2) 

\z-z'\ = 2r sm(n(9 - 9')/2) . (88) 

One thus obtains 



SV e SV e ' 

J y n y n 

v 2 



SA^evtsJoi ^lO'oi 



expf cosn(9 - 9') ) J [ 2j 01 sin 



W 

expf -^L coan (e - 9') ) Jo ( 2 joi cos 



2x 2 



Setting 9' = 9, we obtain the standard deviation on T^f {oo} 
<W) 2 ) _ 1 



V 2 



2N cvtsJ 2 1 j 1 (j 01 y 



r 



n(9 



(89) 



(90) 



These relations show that the relative statistical error 
on integrated flow depends on the parameters V n and a 
only through the resolution parameter x, as explained in 
Sec. VII A. One sees that the error diverges exponentially 
for small x> while the divergence is only polynomial with 
finite-order cumulants [34]. 

The linear correlation coefficient between two esti- 
mates V^f{oo} and V% {oo} is defined as 



SV*5V? 



(91) 



It may vary between — 1 and 1 , which correspond to the 
cases when SV^{oo} and 6V~ {oo} are opposite or equal, 
respectively. The variation of the correlation coefficient 
c{9 -9'), computed with the help of Eqs. (89) and (90), 
is displayed in Fig. 6 as a function of the relative angle 
for several values of x- One sees that the correlation is 
maximum for 9' — 9 and 9' = 9 + n /n, which is natural 
since V® +7T ^ n = V®. For small values of x, the correlation 
vanishes when the relative angle 6' — 9 is large enough, 
while for larger values of %, an anticorrelation appears 
around 9' = 9 + n/2n. 



The shape of the correlation function c(9 — 9') shows 
that different estimates (for different values of 9) are not 
necessarily strongly correlated. In order to decrease the 
statistical error, one can thus average V% {oo} over sev- 
eral values of 9. These can be chosen equally spaced, 




FIG. 6: Correlation function c(0 - 9'), denned by Eq. (91), 
as a function of n(8 — 9') for three values of the resolution 
parameter, % — 0-5, 1, 1.5. 



i.e., 9 = kir/p with k = 0, ■ • ■ ,p — 1. The statistical 
error on this average can be derived from Eq. (89). Nu- 
merical estimates for various values of x an( i different 
numbers of values of 9 are given in Table II, where we 
have assumed that JV ev t s = 20000 events have been ana- 
lyzed. This Table shows two features. First, the value of 
the error badly diverges when the resolution parameter 
X decreases, so that it is very unlikely that our method 
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can be applied when \ is smaller than 0.5, unless a huge 
number of events are available, in agreement with our 
discussion in Sec. VII D. Second, the statistical error de- 
creases when the number of values of 9 increases, but 
quickly saturates, so that in practice, 4 or 5 values of 9 
are enough. 



TABLE II: Relative statistical error on V n {oo}, for N cvtB — 
20000 events and various values of the resolution parameter 
X and the number of values of 8. 



Nb of points 


X = 0.6 


X = 0.7 


X = 0.8 


x=i 


X = 1.5 


1 


22.2% 


7.7% 


3.8% 


1.70% 


0.75% 


2 j 


15.7% 


5.4% 


2.7% 


1.18% 


0.46% 


3 


12.8% 


4.4% 


2.2% 


0.98% 


0.41% 


4 


11.4% 


4.0% 


2.1% 


0.94% 


0.41% 


5 


11.0% 


3.9% 


2.0% 


0.94% 


0.41% 


+oo 


10.9% 


3.9% 


2.0% 


0.94% 


0.41% 



The statistical error is of course minimum when the 
number of values of 9 goes to infinity. It is then given by 
averaging Eq. (89) over both 9 and 9'. One obtains 



({5Vnf) 

V 2 



ioi 



(92) 

Values in the last line of Table II were calculated using 
this expression. 

It is instructive to compare the statistical error on 
V^{co} to that on estimates from cumulants of 2k- 
particle correlations (where k = 1,2, ■■■), denoted by 
V n {2k} in Ref. [34]. In both cases, the statistical er- 
ror depends on two parameters, namely the total num- 
ber of events, N evts , and the parameter x- However, 
the x-dependence of the errors is very different. Thus, 
whereas the statistical uncertainties on the flow estimates 
from cumulants for small x diverges like a power law 
l/x 2fc (see Appendix D in Ref. [34]), the error on l^{oo} 
depends exponentially on \ [ sec Eq- (90)]. That ex- 
plains why the detected multiplicity per event, which 
influences the value of Xi pl&y s a most crucial role in 
the present method, while cumulants could accommodate 
lower statistics. 

A numerical comparison is made in Table III for vari- 
ous values of the resolution parameter x- As anticipated, 
for small Xi the relative error on V^{oo} is significantly 
larger than that on V n {2k}. One also notes that the sta- 
tistical error on V^{oo} is always larger than the error 
on the standard flow estimate (from two-particle corre- 
lations) V"„{2}. For values of x of order unity, however, 
the statistical error on 14,{oo} is not much larger than 
the error on the two-particle estimate V^{2}. One must 
keep in mind here that lower order estimates such as 
I4{2} are biased by nonflow correlations, and that the 
resulting systematic error is usually much larger than the 
statistical error. 



TABLE III: Comparison between different methods: the rel- 
ative statistical error on the integrated flow V„ is shown for 
the cumulant method [34], and various cumulant orders (de- 
noted by V n {2k} with integer k), and for the present method 
(V^{co}). As in the previous table, the number of events is 
Acvts = 20000 events, and the resolution parameter x takes 
several values. 





X = 0.6 


X = 0.7 


X = 0.8 


x=i 


X = 1.5 


5V n {2}/V n 


1.3% 


1.0% 


0.83% 


0.62% 


0.37% 


5V n {4}/V„ 


4.5% 


2.7% 


1.8% 


1.00% 


0.43% 


5V n {6}/V„ 


7.7% 


3.7% 


2.1% 


0.99% 


0.41% 


5V n {8}/V„ 


9.9% 


4.1% 


2.1% 


0.95% 


0.41% 


5V n {oo}/V n 


10.9% 


3.9% 


2.0% 


0.94% 


0.41% 



One also sees clearly in Table III that the error on 
the present estimate T^„{co} is the limit of the error on 
V n {2k} as k goes to infinity, for fixed x- This reflects 
the fact that T4,{co} is itself the limit of V n {2k} as k 
goes to infinity, as explained in Sec. III. It is interesting 
to note that the statistical error does not increase mono- 
tonically with the order of the cumulant 2k, as one might 
think. This was already pointed out in Ref. [34]. For a 
fixed Xi t ne error on V n {2k} first increases as k increases, 
reaches a maximum for some value of k and then slightly 
decreases. This can be seen in the last four columns of 
Table III. While the statistical error is always smallest for 
the standard, two-particle estimate V n {2}, the error on 
Kj{oo} is smaller than that on the estimate from four- 
particle correlations 14,(4} as soon as x exceeds 0.89. 
Comparing with the values of \ given in Table I, this 
suggests that the new method can be used without prob- 
lem at RHIC. The statistical error on our flow estimate 
is expected to be comparable to that on the estimate 
from four-particle cumulants, and significantly smaller 
than the systematic error (due to nonflow correlations) 
on the estimate from the standard method [24] . 



F. Statistical error on the differential flow 

We now turn to the statistical uncertainty on the dif- 
ferential flow estimates v'^ n and its average over 9, v' mn . 
Note that, unlike the case of integrated flow, the errors 
we shall compute in this Section are absolute errors, not 
relative ones. 

Since differential flow is analyzed in principle in a nar- 
row phase-space window, the statistical error on differen- 
tial flow is much larger than the error on integrated flow. 
We therefore neglect the statistical error on integrated 
flow here: in Eq. (12), we replace y„{oo} and r® by their 
theoretical values, i.e., V n and ro [see Eq. (24)]. Simi- 
larly, we replace in the denominator the sample average 
by an expectation value: 



/giroC}' 
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dG cA .(ir) 
dr 1 



(93) 



becomes 



where the rhs is given by Eq. (86). Equation (12) then 



j 2 /(4 X 2 ) 

«£„{«>} = (-ir X(— y { cos ("^ - )) Rc ( imG ( r o } • 



(94) 



We now evaluate the correlation between t^„{oo} and w^„{oo}, using an equation similar to Eq. (70), with N cvts 
replaced by the total number of "differential particles" , N' . We assume, for simplicity, that ip and G(z) are statistically 
independent (which means in particular that the differential flow vanishes) so that the term involving ip factors out. 
Averaging over -0, this term gives 



{cos{mn(ip — 9)) cos(mn(ip — 8'))) = — cos(mn(9 — 6')). 



One thus obtains 



6v% n {co}6vZ n {co}) = 



w> 



cos(mn(9-9'))e 3 «^ ( - 2x 
2N'J m (joi) 2 



Re(i m e iroQl> )Re(i m e iroQe )) . 



(95) 



(96) 



Now, Re(i m G(zj) is, up to a sign, the real (resp. imaginary) part of G(z) for even (resp. odd) m. Using Eq. (68) 
and Eq. (23), one finally obtains 



*«lW 6v mn{°°} 



cos(mn(9 — 9')) 
4iV'J m (ioi) 2 



exp| Ml cos (n(9- 6'))^ J (2j 01 sin( ^ ^ 
+ (-l) m cxp^-^cosn(0-0')) Jb^joi cos^ , ; 



, (97) 



where N' is the number of particles in the consid- 
ered phase-space bin. One checks that this expression 
is invariant under the transformation 9 — ► 9 + n/n, 
as expected from the symmetry property w^„{oo} = 

v'rntS^'ioD}. Setting 9 = 6' in Eq. (97) yields the ab- 
solute statistical error on v'^ n . 

As in the case of integrated flow, one must average 
over several values of in order to reduce the statistical 
error. Tables IV and V show the variation of the error 
as the number of values of 9 increases, for the first two 
harmonics v' n {oo} and f2„{oo}. Once again, 4 or 5 values 
of 9 are enough in practice to minimize the statistical 
error. In the limit when the number of 9 goes to infinity, 
the error is given by 



«m{°°} 2 ) - Wmnf = 

E ^Ooi) 2 W(0). (98) 



2A'J m Ooi) 2 ^ J9(j ° l)2/9+m (5 



It is also interesting to compare the error computed in 
this Section with the statistical errors on differential flow 
estimates from multiparticle cumulants [34]. This com- 
parison is performed in Tables VI and VII for the first two 
harmonics v' n and v' 2n , respectively. The same comments 
apply here as for the integrated flow (Sec. VII E): statis- 



TABLE IV: Statistical error on the reconstructed value of the 
differential flow v' n , for N' = 6 x 10 5 particles in the bin, and 
various values of the resolution parameter and the number of 
values of 9. Note that these are absolute values, not relative 
values as in Tables II and III. 



Nb of points 


X = 0.6 


X = 0.7 


X = 0.8 


x = i 


X =1.5 


1 


6.9% 


2.4% 


1.19% 


0.53% 


0.24% 


2 


4.9% 


1.7% 


0.84% 


0.37% 


0.17% 


3 


4.0% 


1.4% 


0.69% 


0.31% 


0.14% 


4 


3.5% 


1.2% 


0.63% 


0.29% 


0.14% 


5 


3.4% 


1.2% 


0.62% 


0.29% 


0.14% 


+oo 


3.3% 


1.2% 


0.62% 


0.29% 


0.14% 



tical errors on ^^{oo} are much larger than errors on 
standard flow estimates for small x, but become compa- 
rable, and even smaller than estimates from higher-order 
cumulants, when \ is close to unity or larger. 

In addition, one notes that the error on the higher 
harmonic v' 2n is only slightly larger than the error on v' n , 
while in standard flow analyses the error on higher har- 
monics is often much larger. It turns out that in the case 
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TABLE V: Same as Table IV, but for the higher harmonic 

«2n- 



Nb of points 


X = 0.6 


X = 0.7 


X = 0.8 


x = i 


X=1.5 


1 


8.3% 


2.9% 


1.43% 


0.63% 


0.28% 


2 


5.9% 


2.0% 


1.02% 


0.46% 


0.22% 


3 


4.8% 


1.7% 


0.83% 


0.37% 


0.17% 


4 


4.1% 


1.4% 


0.72% 


0.32% 


0.15% 


5 


3.8% 


1.3% 


0.69% 


0.32% 


0.15% 


+oo 


3.7% 


1.3% 


0.68% 


0.32% 


0.15% 



TABLE VI: Comparison between different methods: the sta- 
tistical error on the differential flow v' n is shown for the es- 
timates from cumulant of 2-particle (v' n {2}) and 4-particle 
(«n{4}) correlations [34] and for the present method (v' n {oo}). 
As in the previous table, we assume N' = 6 x 10 s particles in 
the bin, the resolution parameter \ takes several values. 





X = 0.6 


X = 0.7 


X = 0.8 


x = i 


X=1.5 


5v' n {2} 


0.18% 


0.16% 


0.15% 


0.13% 


0.11% 


5v' n {4} 


0.88% 


0.61% 


0.45% 


0.29% 


0.15% 


Sv' n {co} 


3.3% 


1.2% 


0.62% 


0.29% 


0.14% 



of higher harmonics, the statistical error on v' mn {oo} is 
even smaller than the statistical error on the estimate 
from the standard flow analysis [which we denote by 
v' mn {m + 1} since it involves a (m + l)-particle correla- 
tion], if the resolution parameter \ is large enough. More 
precisely, our method yields smaller error bars than the 
standard method as soon as x > 1-3 for v' 2nl %> 1.1 for 
«3„, and x > 0-94 for v' 4n . Even from the point of view 
of statistical error bars, it seems to be the best method 
to measure higher harmonics W3, W4, or at least to set up- 
per bounds on their values, which could easily be done 
at SIS or AGS energies where the reaction plane resolu- 
tion is high: correlating more particles may in some cases 
lead to smaller statistical fluctuations, which is a rather 
unexpected conclusion. However, one must keep in mind 
that systematic errors on these higher harmonics may be 
large, as shown in Sec. VB. 



VIII. DETECTOR EFFECTS 

Until now, we have performed all calculations assum- 
ing that the detector is perfect, in the sense that it is 
azimuthally isotropic. In this Section, we shall relax this 
assumption, and investigate the influence of acceptance 
inefficiencies on the determination of flow through the 
present method. As we shall see, the generating function 
will acquire a phase, which does not affect the fact that 
the zeroes lie on the imaginary axis. The only sizeable ef- 
fect on the ^-averaged estimate U„{oo} is of second order 



TABLE VII: Same as Table VII, but for the higher harmonic 
v' 2n . Errors on estimates from 3- and 5-particle cumulants [34] 
are compared with the present estimate v' 2n {oo}. 





X = 0.6 


X = 0.7 


X = 0.8 


x = i 


X = 1.5 


Sv' 2n {3} 


0.48% 


0.38% 


0.32% 


0.24% 


0.16% 


Sv' 2n {5} 


1.42% 


0.88% 


0.59% 


0.33% 


0.16% 


Sv' 2n {oo} 


3.7% 


1.3% 


0.68% 


0.32% 


0.15% 



in the acceptance coefficients a n which we soon define, 
and is negligible in most practical cases. 

Let us denote by A{4>) the probability that a particle 
with azimuthal angle cj) be detected: A(cf)) represents the 
acceptance-efficiency profile of the detector. We choose 
the normalization A(4>) d4>/2ir = 1. 17 We assume for 
simplicity that A{cf)) is independent of the rapidity and 
transverse momentum of the particle. Since A(4>) is a 
27r-periodic function, it is natural to expand it in Fourier 
series. We denote by a p the corresponding Fourier coef- 
ficients: 

a n = JJe^A^)^-. (99) 

The normalization choice of course reads do = 1. 

The probability distribution of for particles seen in 
the detector, for a fixed orientation of the reaction plane 
is then 

=£=m E Vne in ^-**\ (100) 

^ n— — 00 

where we use the conventions V- n = v n , Vo = 1. 

Using this distribution, we can compute the generating 
function as in Sec. IIIC. Equation (17) still holds, but 
(Q 6 \&r) is no longer given by Eq. (18), which holds only 
with a perfect acceptance. From the definition of Q e , 
Eq. (3), we obtain instead 

+00 

(Q 6 \$r)= E V m Rc(a n - m e mS e -^«), (101) 

m— — 00 

which replaces Eq. (18) for an arbitrary detector. For a 
perfect acceptance (ao = 1, a n ^o = 0), this expression 
reduces to Eq. (18), as should be. 

We shall now assume for simplicity that only one 
Fourier harmonic of the flow, v n , is non- vanishing, with 



Strictly speaking, such a normalization amounts to assuming 
that the overall efficiency is 100%, so that it can be larger than 
100% for some values of the azimuth if the acceptance is not 
isotropic. . . However, normalizing A((j>) a priori does not influ- 
ence the discussion. This simply reflects the fact that the flow 
analysis is not sensitive to the overall efficiency of the detector, 
but only to the anisotropics in the acceptance. 
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n 7^ 0. This is most probably a reasonable assumption 
at ultrarclativistic energies. Then, only three terms con- 
tribute in the sum, m = —n, m = and m — n. We may 
then rewrite the previous equation as 

(Q 8 \$r) = V Re(a n e in6 ) 



+ V n Re 



(1 + a 2n e 2ln9 ) e m ^ R - 0) , (102) 



where Vo = (^2 w j) an d we have used the normalization 
condition do = 1, as well as the properties a^ 2n = a 2r n 
see Eq. (99), and V- n = V n , which follows from Eqs. (2) 
and (4). Inserting this value in Eq. (17), we obtain 



e zQ ° |$«) = exp ( z V Re (a n e in6 ) + + z V n Re f (l + a 2n e 2md ) e 1 ^^ 



(103) 



Averaging this expression over (we still assume that 
a is independent of we obtain 

G\z) = e zVaK < a " eme ^ 2 * 2 / 4 I (zV n |1 + a 2n e 2me \) . 

(104) 

Comparing this result with the perfect-acceptance case 
in Eq. (22), there appear several modifications. First of 
all, G (z) is no longer real- valued for values of z on the 
imaginary axis. Then, there is a new term in the expo- 
nential factor, which involves the acceptance coefficient 
a n . This term cannot produce any spurious zero of G e (z) , 
nor shift the position of the zeroes. When z is purely 
imaginary, it is a pure phase, which goes away when tak- 
ing the modulus |G e (,2)|. However, this phase makes the 
real part of G e {z) oscillate, and thus makes it vanish on 
the imaginary axis. This explains why it was important 
to consider (the minimum of) the modulus, not (the zero 
of) the real part, as mentioned in Sec. IIIC. 

The remaining modification is more relevant for the 
procedure introduced in Sec. II B, since it results in a 
shifting of the position of the generating-function zeroes. 
The argument of the Bessel function in Eq. (104) now 
depends on 9, and involves the acceptance coefficient a 2n . 
Thus, using Eq. (9), we now obtain 



V°{oc} =V n \l + a 2n e 



2in8 I 



(105) 



where a 2n is given by Eq. (99): the estimate V% now 
depends on 0, as can be seen in Fig. 4. As expected from 
the general discussion in Sec. II B, it is (7r/n)-periodic. In 
the case of a detector with a hole of a radians around <fi = 
it, an explicit calculation gives a 2n = — sin(na)/[(27r — 
a)n]. Thus, in the simulation of Sec. II D, ~ 0.0827, 
quite a small number, but large enough to yield the 9- 
dcpcndcnce observed in Fig. 4. 



Averaging Eq. (105) over 9, one obtains 



KM = K(l + |a 2 „| 2 ). 



(106) 



The relative error is 



a 2r , 



which is in practice small if 



the detector has a reasonable azimuthal coverage. 

Consider now differential flow. We leave open the pos- 
sibility that the corresponding acceptance function is not 
the same as for integrated flow, and we denote by a' n the 
Fourier acceptance coefficients for the differential parti- 
cle. For the sake of brevity, we shall only consider the 
case m — 1 (differential flow in the same harmonic as in- 
tegrated flow); the generalization to m arbitrary only in- 
volves much more tedious calculations. Finally, we make 
the same assumption as above, namely, that only one 
flow harmonic does not vanish. 

Under this assumption, Eq. (36) is replaced with an 
equation similar to Eq. (102) 



(cos(n(^ - = v' n Re [(1 + a' 2n e 2me )e m ^-^ 

(107) 

This identity can be then combined with Eq. (103) to 
compute D B m {z) defined by Eq. (30). To perform the 
integration over one can use the following integral: 



2 \c (Be im *) e z K < Ae ^ ^ = MB A f - l,„ ( : | . l| ). 



2tt 



\A\r 



(108) 

valid for arbitrary complex numbers A and B. We use 
this result with m = 1 (other values of m might be useful 
for higher harmonics) and <p = $_r — 9: 



!>!(*) -e |l + a 2n e M "*| h{zV n \1 + a 2n e \) 



(109) 



r 



Note that the exponential prefactor and the argument of the Bessel function are the same as in Eq. (104). After 
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some algebra, similar to that in Sec. HIE, one finally 
obtains 



C{oo} 



, Re[(l + 4 n e 2 ^)(l+aLe- 2 ^)] 
|l + a 2 „e 2 ^| 



• (no) 



As in the case of integrated flow, acceptance inefficien- 
cies result in the appearance of a ^-dependent, (ir/n)- 
periodic, multiplicative coefficient in front of v' n . In the 
particular case where a' 2n = a 2n , Eq. (110) gives 

v'°{™} = v' n \l + a 2n e 2me \. (Ill) 

Integrating over phase space, one recovers Eq. (105). 

IX. CONCLUDING REMARKS 

In this paper, we have introduced a new method of 
analysis of anisotropic flow in heavy-ion collisions. This 
is the first method to extract both integrated and dif- 
ferential flow from the interparticle correlations between 
a large number of particles, rather than from the corre- 
lations between only a small number. Thus, it is able 
to isolate a collective effect which involves all particles, 
such as flow, removing the influence of other, "nonflow" 
effects, which can bias few-particle methods. The flow 
estimates obtained with the new method are therefore 
more reliable than with any other method. 

This increased reliability is of course important inas- 
much as flow is concerned: with more accurate col- 
lective flow estimates, the constraints on models be- 
come stronger. But it also matters for the measure- 
ment of other effects, as for instance when reconstruct- 
ing jets from the azimuthal correlations they induce be- 
tween high momentum particles at ultrarelativistic ener- 
gies [16, 25, 54]. At high pr, (elliptic) flow is large, and 
constitutes a huge background. It is most important to 
know its magnitude accurately to extract nonflow corre- 
lations. 

Let us enter into a little more detail. The method 
relies on the search of the first minimum of a generat- 
ing function of multiparticle azimuthal correlations; the 
actual position of this minimum directly yields an inte- 
grated flow estimate. Once the minimum has been found, 
computing a second function at its position gives differ- 
ential flow. This procedure is simpler to implement than 
other methods of flow analysis. It is significantly faster 
numerically, and less tedious, than the methods based 
on cumulants of multiparticle correlations. In addition, 
in the differential flow analysis there is no need to sub- 
tract autocorrelations, which are negligible. Finally, de- 
tector effects are negligible as well in most cases. We have 
carefully studied the various sources of errors, either sys- 
tematic errors intrinsic to the method (due to nonflow 
effects, higher flow harmonics, detector inefficiencies), or 
statistical errors. This allows us to conclude that it is 
the most accurate method to measure both directed and 
elliptic flow in collisions from 100 MeV to a few GeV per 



nucleon (SIS and AGS energies) and elliptic flow at ultra- 
relativistic energies, especially at RHIC and the CERN 
Large Hadron Collider (LHC). 

It may be possible to extend the method, so as to be 
able to measure directed flow at ultrarelativistic energies 
as well, paralleling the step which led from the cumulant 
method of Rcfs. [33, 34] to that of Ref. [35]. Namely, the 
generating function of one complex variable G{z) could 
be generalized to a two-variable function G(zi,z 2 ) in- 
volving azimuthal correlations in two different harmonics 
(corresponding to V\ and v 2 ) at once. This new function 
could then be used to measure directed flow using as a 
reference the elliptic flow, which is large at ultrarelativis- 
tic energies, and thus provides a good reference. 

Another generalization of the method would be to ex- 
tract directly nonflow azimuthal correlations, and es- 
pecially, to obtain them with respect to the impact- 
parameter direction in non-central collisions. This could 
be done, for instance, to study the azimuthal depen- 
dence of Hanbury-Brown Twiss (HBT) correlations [55], 
or that of jet quenching [54]. 

Finally, thanks to the generality of its formalism, this 
method could easily be extended to other types of ob- 
servables where "fluctuations" are of interest. It seems 
to be the most natural method to look for critical, large 
scale fluctuations in a system, which are expected in the 
vicinity of a phase transition [47] . 
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APPENDIX A: AN ALTERNATIVE FORM OF 
THE GENERATING FUNCTION 

The generating function on which our method is based, 
G e (z), has been defined in Eqs. (3) and (6): 

G e {z) = |exp^5Z«;jCos(n(^ -fl)) | } . (Al) 

An alternative choice is: 

{M 
JJll + ^cosK^ -«))]>. (A2) 
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The procedure to determine the integrated flow with 
G e (z) is exactly the same as for G e (z), as presented in 
Sec. II B. In order to determine the differential flow, 
one should replace e 4r ° < 3 in the numerator of the rhs 
of Eq. (12) with the term in curly brackets in the rhs 
of Eq. (A2), with z = ir®; next, one should replace 

in the denominator of the rhs of Eq. (12) 

with the derivative dG s /dz, evaluated at z = ir$. 

The essential point is that & (z) factorizes for indepen- 
dent subsystems, as does G e (z), so that the discussion in 
Sees. IIIB and IIIC can be readily extended to G e (z). 
The only difference is that a in Eq. (17) is no longer given 
by Eq. (19). In fact, a vanishes for independent particles, 
and is generally smaller for & than for G e . Although the 
value of the generating function is modified, the zeroes 
remain at the same position, and in particular the first 
one, from which flow (cither integrated or differential) 
is determined. A similar equivalence is also known in 
statistical mechanics [56] . 

A generating function similar to G (z) was chosen in 
the cumulant expansion of Ref . [34] . The argument was 
that when expanding in powers of z, G e (z) involves only 
correlations between different particles, while G e (z) in- 
volves "autocorrelation" terms. The bias induced by 
these autocorrelations was studied in detail in Ref. [33] 
for finite-order cumulants. In general, this bias is of the 
same order of magnitude as the bias induced by nonflow 
correlations, as explained in Sec. Ill E 4. 

Let us illustrate this by a simple explicit example. We 
repeat the calculations of Sec. IV with G e {z) instead of 
G®{z): particles are emitted in collinear jets containing q 
particles each. If the jet angles are randomly distributed, 
one obtains 



G d (z) = ((l + 



zcost 



K M/q 



(A3) 



where angular brackets denote an average over <ft, which 
gives: 



[9/2] 



.,21 



M/q 



^ (g-2Z)!(Z!) 2 2^ 



(A4) 



1=0 



For independent particles (q = 1) G e (z) — 1 [compare 
with Eq. (44)]. This function has no zero, hence the 
analysis yields no spurious flow. This shows that G e (z) 
is an improvement over G s (z) when there is no nonflow 
correlation. For higher values of q, G e (z) is a polynomial 
whose roots can be computed numerically. The estimate 
V^{oo} is then given by Eq. (9). One obtains the values 
1.70, 2.94, 4.07 for q = 2,3,4 respectively. This is very 
close to the value y„{oo} = q obtained with G e (z) [see 
Eq. (46)], which shows that G (z) no longer represents 
an improvement over G e (z) when nonflow correlations 
are present. 

The same conclusions hold for systematic errors, whose 
order of magnitude is not modified by autocorrelations. 



In practice, the latter may increase the contribution of 
each term in the rhs of Eq. (48), but this equation re- 
mains valid as an order of magnitude, with cither form 
of the generating function. 

The final results of Sees. V to VIII also apply when 
G 6 {z) is used instead of G e (z), although intermediate 
calculations may differ. In particular, a derivation of 
statistical errors with a generating function similar to 
G e (z) is given in Appendix D of [34]. 

As a conclusion, G e (z) is expected to give more ac- 
curate results than G e (z) when nonflow correlations are 
weak. The price to pay is that G e (z) requires much more 
computer time: for each value of z, one needs to evalu- 
ate the product over all particles in Eq. (A2), instead of 
calculating only one flow vector, valid for all values of z. 
Furthermore, since the motivation of this paper was to 
obtain a method which remains valid when nonflow cor- 
relations are present, and since G (z) is not better than 
G 6 {z) in this respect, we chose to use G (z) in the paper. 



APPENDIX B: RELATION WITH THE 
CUMULANT EXPANSION OF REF. [33] 

In Ref. [33], integrated flow was determined from the 
following generating function (Eq. (B2) of Ref. [33]): 



G( Zl ,z 2 ) = {e^ Q+Z2Q '}, 



(Bl) 



where z\ and zi are two independent complex variables, 



Q = Q x + iQ y , and Q* 



Cumulants 



((Q k {Q*) 1 )) were then defined by the power-series ex- 
pansion of lnG(zi, z 2 ): 



Z k zh 



k.l 



(B2) 



The new generating function, Eq. (6), can be expressed 
in terms of G(z\, z 2 ): 



(B3) 



The cumulant (((Q e ) k )) defined by Eq. (15) is a linear 
combination of the ((Q fe "'(Q*)')): 



««m> = iE 



2 k ^ \ I 

i=o x 




W-W ((Q k - l (Q*) 1 )). (B4) 



In Ref. [33], the recommended procedure was to extract 
an estimate of the flow, denoted by V n {2k} in Ref. [34], 
using the "diagonal" cumulants ((Q k (Q*) k )), which are 
the only non-vanishing cumulants except for statistical 
fluctuations and detector asymmetries. By definition of 
the estimates, ((Q k (Q*) k )) is equal to (V n {2k}) 2k , up to 
a multiplicative factor. 

In the present paper, we obtain an estimate V n from 
{((Q e ) 2k )), which was denoted by V°{2k} in Sec. IIIC: 
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By definition, «(Q e ) 2fe )) is equal to {V°{2k}) 2k , up to 
a multiplicative factor. In order to obtain the relation 
between V*{2h} and V n {2k}, let us average Eq. (B4) 
over 9: 



Jo 



dO 



One then obtains the following relation: 



d9_ 

2tt 



W{2k}y* = (V n {2k}) 



2k 



(B6) 



In the limit of an infinite number of events, and with a 
perfectly symmetric detector V®{2k} is independent of 9 
by azimuthal symmetry. The above equation shows that 
V® {2k} then coincides exactly with the estimate V n {2k}. 
In this case, our estimate V^{oo} is exactly the large- 
order limit of V n {2k} for large k. 

However, this is no longer the case when statistical 
fluctuations and detector asymmetries are taken into ac- 
count. In this paper, we first let the order k go to infinity 
for a fixed value of 9, and obtain the limit V^f{oo}. Then, 
we average over 9. On the other hand, in the method of 
Ref. [33], one first averages over 9 for a fixed k, which is 
natural when working at a finite cumulant order. 

This is actually a minor difference. The only prac- 
tical consequence is that acceptance corrections, which 
must be applied when the detector has limited azimuthal 
coverage, differ. The present approach, where the whole 
analysis is done for a fixed value of 9, is more natural and 
leads to simpler acceptance corrections (see Sec. VIII) 
than the earlier method (for which the acceptance cor- 
rections were derived in Appendix C of Ref. [34]). 

Another difference with the method of [33] is that, in 
the latter, weights Wj in Eq. (2) were chosen propor- 
tional to 1/ VM, while they are independent of M in the 
present paper. The rationale was that lower-order cu- 
mulants contain large trivial contributions from "auto- 
correlations" : for instance, (|Q 2 |) is nonvanishing even 
if there is no flow. With the chosen weight, these con- 
tributions were independent of M and could easily be 
subtracted. However, the magnitude of these autocor- 
relations becomes smaller and smaller as the cumulant 
order increases, and they are negligible with the present 
method (which is the limit of asymptotically large order), 
so that \j\J~M factors are no longer required. 

In Ref. [34], an improvement over the method of 
Ref. [33] was proposed, using a generating function sim- 
ilar to G e (z) in Eq. (A2). This new formulation was 
free from autocorrelations. However, we pointed out that 
when the detector has limited azimuthal coverage, mul- 
tiplicity fluctuations induce errors on lower-order cumu- 
lants which are similar to the errors induced by autocor- 
relations. These errors were minimized by taking weights 
proportional to 1/M, and using a slightly modified defi- 
nition of the generating function of cumulants (Eq. (7) of 
Ref. [34]). With the present method, such a refinement 
is not necessary. 



In this paper, we choose to work with weights which are 
independent of M for two reasons: first, this simplifies 
the discussion of Sec. Ill and makes the relation to Lee- 
Yang theory clearer. Second, the multiplicity may not 
be a relevant quantity in heavy-ion experiments at lower 
energies, where a nuclear fragment contributes as much 
to the flow as several nucleons. However, our analysis 
could as well be performed with M-dependent weights. 
The only difference is that the integrated flow V n would 
scale differently with the multiplicity, and its fluctuations 
with impact parameter, studied in Sec. VI, would also 
differ. 



APPENDIX C: LARGE ORDERS 

Consider a function f(z) which is analytic in the vicin- 
ity of the origin. It can be expanded in power series: 



f(z) = J2a k z k . 



(CI) 



fe=0 



The coefficients a k can be expressed in integral form: 

f(z) dz 



ak = 



z k+1 2?7r' 



(C2) 



where the integration contour circles the origin in the 
complex plane (full line in Fig. 7). The integration con- 
tour can then be expanded until it reaches a singularity 
of f(z) (sec Fig. 7). We are interested in the asymp- 
totic behavior of the coefficients a k as k goes to infinity. 
For large k, the integral in Eq. (C2) is dominated by the 
smallest values of \z\. Therefore, the large-order behavior 
of dfc is determined by the singularities of f(z) which are 
closest to the origin. 

The first case of interest in this paper (Sec. IIIB) is 
f(z) = lnG 8 (z), where G e (z) is analytic in the whole 
complex plane (entire), and satisfies the symmetry rela- 
tion Eq. (8) (i.e., it is real for real z). Let us denote by 
zq the zero of G(z) which is closest to the origin in the 
upper plane. We assume that it is a simple zero. Thanks 
to Eq. (8), its complex conjugate Zq is also a zero. The 
singularity of f(z) at these points is logarithmic, and the 
discontinuity of f(z) across the cut (Fig. 7, left) is — 2iir. 
The integral in Eq. (C2) then reduces to 



ak 



= --Re 



dz 

z k+l 

2_ / 1 



dz 



k \{z ) k J' 
It is real, as expected since G (z) is real for real z. 



(C3) 



The second case we are interested in (see Sec. HIE) 
is f(z) = D 6 m (z) I G e (z) , where G (z) is the same func- 
tion as above, and D 6 m {z) is another function sharing the 
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FIG. 7: Integration contours in the case where f(z) has a 
logarithmic singularity (left) or a pole singularity (right). The 
full line is the initial integration contour, and the dashed line 
is the deformed contour. 



same properties: D e m {z) is entire and also satisfies the 
symmetry relation Eq. (32). The singularities of f(z) 
which are closest to the origin are Zq and z$, as in the 
previous case. They are now pole singularities, and their 
contribution to the integral in Eq. (C2) is given by the 
theorem of residues: 



D e m (zo) 



= -Re 



(z ) fe+ W(*o) (* *) fe+ W(*o*) 
2 D s m (z ) 



(*o) fe+1 (G )'(z o )J ' 



(C4) 



where (G 6 )' denotes the derivative of G e (z) with respect 
to z. 



APPENDIX D: A QUANTITATIVE ESTIMATE 
OF SYSTEMATIC ERRORS 



Because of flow itself, and even in the case of a perfect 
detector, the standard deviation a defined in Eq. (19) 
depends on the orientation of the reaction plane In 
this Appendix, we shall explicitly compute the system- 
atic error on integrated and differential flows due to this 
dependence, assuming that particles are independent (no 
nonflow correlations). In addition, we shall show that 
this dependence is influenced by the higher flow harmonic 
«2n, which thus contaminates the measurement of v n . 
Throughout the Appendix, we assume that the detector 
has perfect azimuthal isotropy. 

For a fixed orientation of the reaction plane the 
normalized azimuthal distribution is 



dN 1 V s r (a 



$7 



(Dl) 



where we define vq = 1 and v_„ = v n . 

In order to compute the generating function (6) in var- 
ious physical situations, it is useful to first evaluate the 
following single-particle average with the azimuthal dis- 
tribution (Dl): 



'dN 



} z cos(n((p—9)) 



1 

+ oo 



= v m n I m {z) cos(mn($ fl - 6)). 

m— — oo 

(D2) 

Please note that both m and — m contribute equally to 
the sum. 

Let us now consider a simple model of the collision: 
we assume that the azimuthal angles of the particles are 
independent for a fixed orientation of the reaction plane 
<f>^, i.e., we assume that all azimuthal correlations arc 
due to flow. If Q e is defined with unit weights in Eq. (3), 
then 



e zQt> \$ R ) = / e zcos ("( < A- e ))|<l), 



M 



(D3) 



Taking the logarithm of this expression, using Eq. (D2) 
and expanding in powers of z up to order z 2 , we obtain 



ln(e zQ I*. 



Mv n cos(n($fl - 6)) z 
- M(v 2n -v 2 n )sin 2 (n(<P R -9))^ 

~2 



+ M (1 + v 2n - 2v 2 n ) — . 



(D4) 



Comparing this result with Eq. (17), one sees that the 
standard deviation a depends on <J>^. This dependence 
was neglected in Sec. Ill C. We shall now assume that it is 
a small correction, and evaluate the resulting systematic 
error on our flow estimates. 

For that purpose, we introduce the following notations: 



Z = Mv n z 
a = n($ R -6) 

= V2n - vj 

e ~ 2Mv 2 n ' 
With these notations, we rewrite Eq. (D4) as 



(D5) 



e M(l+v 2n -2v^)z 2 1 A £ Z cosa-eZ 2 sin 2 a (D6) 



In order to compute the various generating functions, we 
shall need to evaluate the following integral: 



/*27T ] />27T 7 

/ cos(ma)e Zcosa Z 2 sm 2 a^ = - cos(ma) e Zcosa (m 2 - Zcosa) (D7) 
Jo 27T J Q 2n 
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which is obtained with two successive integrations by parts. Using this identity, one derives the following result, which 
is valid to first order in the small parameter e: 



L 



27r 9 9 r/n- f' 2n 

cos(ma)e Zcosa -' z sin a ^ = (1 + em 2 ) / cos(ma) e 

2lT Jo 

= (l + em 2 )I m (Z(l-e)). 



Z(l— e) cos a 



da 
2tt 



(D8) 



r 



Using Eqs. (D6) and (D8) with m = 0, we obtain 



G e (z) = I 
Jo 



= e M(l+V2n-2v n )z /4 /o ( MUnz(1 _ ( D9 ) 

The first zero of G e (z) on the upper half of the imaginary 
axis thus lies at 



ir n = 



^3ol 



Mv n {\ - e) 



(D10) 



The corresponding estimate of the integrated flow is then 
defined by Eq. (9): 



V n {oo}=V n (l-e), 



(Dll) 



where V n = Mv n is the exact result. The relative error 
on the integrated flow is — e. Its expression, Eq. (D5), 
is in agreement with the expected order of magnitude, 
Eq. (48). [Remember that the second term in that equa- 
tion, in l/(Mv n ) 2 , is due to the term in z 3 in the expan- 
sion of ln(e zQ |3>ij). It is thus normal that it does not 



appear in Eq. (Dll), which was derived by truncating 
the power-series expansion to order z 2 .] 

We now evaluate the systematic error on the differen- 
tial flow. Our estimate of v' mn is defined by Eq. (12). We 
assume for simplicity that the particle with angle ip is 
not involved in the definition of Q 6 , i.e., we assume that 
autocorrelations have been removed. 

Since we neglect statistical errors, we can replace sam- 
ple averages by true averages (expectation values) in this 
formula. The denominator of the last factor in Eq. (12) 
is 



dG e {z) ^ 
dz ' z=z ° 

iMv n (l e) e M(i+, 2 „-2^)(4) 2 /4 Ji(joi) _ 

(D12) 



The numerator is evaluated by first taking the average 
value for a fixed which is given by Eqs. (36) and 
(D6): 



^cos[mn(^ - 6)} e z °Q \$ R ) = v' mn e M(i+^„-2<)(^/4 cos(ma) e Z COSa -eZ' .in' a 
where Z = Mv n z®. The average over <& R follows from Eq. (D8) with the help of Eq. (D10): 
/ cos[mn(V - 6)] e^°) - r m v' mn ^(1+^-2^(^/4 ^ + ^ 



(D13) 



(D14) 



r 



Using Eqs. (12), (Dll), (D12), and (D14), we finally ob- 
tain 



v'L{cc} = v' mn (l + em 2 ). 



(D15) 



Therefore, the relative systematic error on v' mn is em 2 , 
where e is defined in Eq. (D5). Note that the correction to 
the integrated flow, Eq. (Dll), has the opposite sign. In 
particular, integrating v'®{oo} over phase space, one does 



not recover the integrated flow. This is because we have 
assumed that autocorrelations have been removed. Note 
also that the correction increases for higher harmonics. 

Finally, note that Eqs. (Dll) and (D15) are quite gen- 
eral: they still hold if there are weights, or if nonflow 
correlations are present. The only difference is that e is 
no longer given by Eq. (D5) (this equation remains valid 
as an order of magnitude, though). 
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